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Linear  Chebyshev  Complex 
Function  Approximation 


I.  Introduction 

The  approximation  of  desired  or  given  functional  behavior  by  finite  sets  of 
simpler  or  specified  basis  functions  is  a  recurrent  problem  in  many  fields.  For 
example,  in  the  mathematical  field,  we  might  wish  to  approximate  a  (desired) 
complex  integral  by  a  set  of  (simpler)  sinusoidal  components.  Or  in  an  antenna 
array  processing  application,  we  often  want  to  realize  a  (given)  low  sidelobe 
behavior  by  means  of  an  array  with  (specified)  element  locations  which  are  not 
under  our  control. 

For  the  case  where  the  given  functional  behavior  and  the  specified  basis  functions 
are  all  real  valued  and  defined  on  a  finite  discrete  data  set,  and  where  the  ap¬ 
proximation  is  afforded  by  a  real-weighted  linear  combination  of  these  basis 
functions,  the  optimum  solution  for  minimizing  the  maximum  magnitude  error, 
i.e.,  the  Chebyshev  norm,  is  in  very  good  shape  due  to  a  fine  algorithm  given  in  [1]. 
Specifically,  this  algorithm  solves  the  following  mathematical  problem:  given  real 
constants  {f;},  { hjk } ,  where  Ki<m,  Kk<n,  m^n,  the  real  quantities  {a,.}^  are 
determined  that  minimize  the  maximum  absolute  value  of  the  error  residuals 

e,  =  F  -  Z  akhjk  for  Ki«m.  (1.1) 

k  =  1 

This  algorithm  has  recently  been  used  to  good  advantage  in  an  array  processing 
application  to  design  some  real  symmetric  weighting  functions  with  very  good 
sidelobe  behavior,  subject  to  constraints  on  the  rate  of  decay  of  the  distant  sidelobes 
[2]. 

Here  we  wish  to  employ  the  algorithm,  as  described  above  for  real  variables  in 
( 1 . 1 ),  for  the  minimization  of  the  Chebyshev  norm  of 

en(z)*f(z)-i  Wz)  (,-2) 

k  -  1 

when  f(z)  and  {hk(z)}y  are  complex,  and  z  can  take  values  in  an  arbitrary  finite 
discrete  point  set.  The  weighting  coefficients  { ak }y,  may  be  complex,  or  alter¬ 
natively,  they  may  be  restricted  to  be  real.  Applications  are  afforded  by  an  antenna 
array  with  arbitrarily  specified  element  locations,  but  employing  weights  that  are 
restricted  to  be  real,  or  alternatively  by  array  weights  that  are  also  allowed  to  be 
phased  (complex).  In  Section  II,  the  basic  mathematical  theory  and  algorithm  for 
the  minimization  of  (1.2)  is  developed.  Numerical  examples  and  applications  of  the 
technique,  some  efforts  attempted  for  extending  the  method  to  a  continuum  of 
values  of  z,  and  a  discussion  constitute  the  rest  of  the  main  body  of  the  report.  An 
Appendix  presents  a  computer  program  in  a  form  w  hich  should  be  useful  to  readers 
interested  in  applying  the  technique  to  their  own  particular  applications. 
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Although  the  above  algorithm  [I]  is  limited  to  a  discrete  set  of  points,  it  has  been 
used  fruitfully  to  minimize  the  continuous  error  (1.2)  over  a  real  variable  z  in  the 
interval  [za,  zb],  when  f  and  {hk}  are  real,  in  the  following  manner.  First,  an  initial 
set  of  m>n  real  points  {z'j1}"1  was  specified  and  the  Chebyshev  norm  minimized  in 
the  usual  fashion,  resulting  in  the  coefficient  set  {a^*}"  For  this  set  of  optimum 
coefficients,  the  locations  of  the  largest  peaks  of  |en(z)|  were  located,  by 

setting  the  derivative  e^(z)  to  zero  and  solving  numerically  for  {z'^}J;  the  number  t 
of  such  peaks  will  generally  be  less  than  m,  but  larger  than  n.  (This  approach 
presumes  the  availability  of  computable  expressions  for  f'(z)  and  {hk(z)}").  Then 
the  modified  set  of  points  { z*?* } J  were  used  for  another  Chebyshev  minimization, 
resulting  in  coefficient  set  {a*“'}".  Repetition  of  this  procedure  stabilized  after  a  few 
trials  with  a  unique  set  of  {z,}{  at  which  the  maximum  errors  were  equal  and 
irreducible.  In  the  examples  tried  in  (2],  the  number  of  peaks,  I,  at  which  the 
magnitude  error  |en(z)|  was  largest  and  equal,  turned  ouf  to  be  n+i.  Further 
discussion  of  this  recursive  approach  is  given  in  Section  V. 

II.  Mathematical  Theory  and  Algorithm 

Let  f  and  hp  .  .  .  ,  hn  be  complex  valued  functions  defined  on  the  finite  discrete 
point  set  Qm  =  {z,,  ....  zm}.  For  a  complex  vector  a  =  (a,,  .  .  .  ,  an)  tCn,  define 
the  complex  error 

^(z>  “  2*  ( akhk(z)  =  en(z;a),  zzQm.  (2.1) 

The  discrete  linear  Chebyshev  approximation  problem  is  to  find  a  complex*  vector 
a  =  (a, . an)  cCn  so  that 

En(f)  =  min  max  |en(z;a)|  =  max  |  en(z;a)  |  .  (2.2) 

atCn  zrQ  ztQm 

m  m 

The  quantity  En(f)  is  called  the  discrete  Chebyshev,  or  minimax,  error  of  the  ap¬ 
proximation  on  the  point  set  Qm. 

We  do  not  solve  this  problem  exactly.  An  algorithm  presented  in  [3]  for  its 
solution  is  erroneous;  we  have  discovered  examples  (see  Section  IV)  such  that  the 
recursive  procedure  described  there  need  not  converge  to  a  solution  of  (2.2).  We  will 
show  that  problem  (2.2)  can  be  replaced  by  a  related  approximate  problem  solvable 
by  available  linear  programming  techniques.  The  exact  solution  of  this  related 
problem  yields  approximate  solutions  of  (2.2).  The  error  in  these  approximate 
solutions  to  (2.2)  can  be  determined  and,  in  fact,  made  arbitrarily  small,  using  the 
results  we  prove  below;  see  Theorems  1  and  2. 

It  can  be  shown  by  standard  mathematical  methods  [4,  p.IJ  that  a  vector  a 
satisfying  (2.2)  exists,  although  it  may  not  be  unique.  Sufficient  conditions  are 
known  that  result  in  unique  a;,  but  we  do  not  need  these  conditions  here.  Therefore, 
no  further  assumptions  on  f,  hp  .  .  .  ,  hn  or  the  point  set  Qm  are  made.  In  order  to 
proceed, we  need  the  following  result. 

*  The  restriction  of  £  to  real  values  is  discussed  at  the  end  of  this  section. 
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Lemma  I.  If  z  =  x  +  iy,  where  x  and  y  are  real,  then 

|z|  =  max  (x  cos  0  +  y  sin  0)  .  (2.3) 

-n<0<n 

Proof.  If  z  =  0,  the  result  is  obvious.  Suppose,  then,  that  z*0.  By  the  Cauchy- 
Schwartz  inequality,  for  every  real  0, 

x  cos  0  +  y  sin  0<  (x2  +  y2)l/2  (cos2  0  +  sin20)l/2  =  |zj 

so  that 


max  (x  cos  0  +  y  sin  0)  <  |z|  . 

-rKS^n 


For  the  particular  value  0  =  arg(z),  it  is  seen  that  (2.3)  holds.  This  completes  the 
proof. 

Now,  let  the  real  and  imaginary  parts  of  the  complex  error  en  (z;a)  be  denoted  by 
Rn(z;a)  and  In  (z;a),  respectively.  Thus,  from  Lemma  1, 

|en(z;a)|  =  max  (Rn(z;a)cos0  +  In(z;a)sin0)  .  (2.4) 

-ir«Kir 

If,  in  this  last  equation,  we  take  the  maximum  over  any  finite  subset  T  of  angles  0  in 
the  interval  -tK0<it,  instead  of  all  angles  in  the  interval  -n<0<n,  we  must  have 

|en(z;a)|  >max  (Rn(z;a)cos0  +  In(z;a)sin0)  .  (2.5) 

0eT 

It  will  be  seen  shortly  that  the  next  result  is  very  important  and  central  to  our 
problem. 


Lemma  2.  Let  6j  =  rr(j- 1  )/p,  j  =  1,2 . 2p,  where  the  integer  p^2.  Let  z  =  x  +  iy, 

and  let 


M  =  max  (x  cos  0.  +  y  sin  0.)  . 
i  =  i . :p  J  1 


(2.6) 


Then 


M  <  |z]  <  M  sec  (n/2p)  .  (2.7) 

Proof.  That  jz|  >  M  is  obvious,  so  we  only  have  to  prove  |z|  <  M  sec  (n/2p).  Let 
P(x,y)  be  the  point  in  the  Euclidean  plane  corresponding  to  the  complex  number 
z  =  x  +  iy  #  0,  so  that 


x  =  |z|  cos  (arg  z) 
y  =  [zj  sin  (arg  z). 
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Thus,  for  any  angle  cp,  we  must  have 

( |z|  cos  (arg  z)  -  x)  cos  cp  +  (  |z|  sin  (arg  z)  -  y)  sin  cp  =  0 

which,  after  simple  algebraic  manipulation,  can  be  written 

|z|  =  x'(cp)  sec  a(cp)  (2.8) 

where  x'(cp)  =  x  cos  <p  +  y  sin  cp  and  a(<p)  =  arg  (ze"1<p).  Alternatively,  (2.8)  can  be 
derived  geometrically  by  considering  x'  to  be  the  x  coordinate  of  the  point  P(x,y) 
after  a  rotation  of  the  axes  through  the  angle  cp.  From  (2.8)  we  have 

|z|  =  x'(0j)  sec  a(0p  ,  j  =  1 . 2p  .  (2.9) 

Let  the  index  k  be  such  that 


M  =  x'(0, )  =  max  x'(0 ). 

14j^2p  J 


(2.10) 


With  the  particular  angles  0^  chosen  here,  x'(0j  +  p)  =  -x'(6j)  for  j  =  1,  .  .  .  ,  p,  so  that 
we  must  have  x'(0k)>O.  Since  z*0  is  fixed  in  (2.9),  it  is  clear  from  (2.10)  and  the 
definition  of  the  angles  <*(0^  that 


Therefore, 


O<seca(0.)  =  min  |seca(0.)|  <  sec(n/2p). 
Kj«2p  J 


I  z  |  =  X'(0k)  sec  a(0k) 

=  Msec(n/2p). 

This  concludes  the  proof. 

We  are  now  in  a  position  to  describe  a  problem  that  we  can  solve  exactly  and  that 
is  related  to  the  given  discrete  linear  Chebyshev  approximation  problem  (2.2).  For 
notational  convenience,  we  define,  for  any  complex  vector  a£Cn, 


Gjfzja)  =  Rn(z;a)  cos  0}  +  In(z;a)  sin  0},  j  =  1 ,  .  .  .  ,  2p  ,  (2.11) 

where  0,,  .  .  .  ,  02p  are  the  angles  given  explicitly  in  Lemma  2.  We  seek  a  complex 
vector  a  =  (a,,  .  .  .  ,  an)  eCn  satisfying 


min 

max 

max  G  (z;a) 

=  1 . 2p  ’ 

a  eC" 

«Qm  J 

max 

max 

G^zid)  . 

(2.12) 

Z£Qm 

j  =  i . 2p 
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Using  standard  mathematical  methods,  it  is  easy  to  see  that  at  least  one  such  vector 
iUCn  exists.  The  connection  between  the  problem  (2.12)  and  the  problem  (2.2)  is 
explored  in  the  next  few  results. 


Theorem  I.  Let  p>2  be  an  integer,  and  let  0j  =  Tr(j-l)/p,  j  =  1,2,  ....  2p.  Then 


M  (f)<  E  (f)<  M 

np'  n'  ' 

np(0  sec  (n/2p) . 

(2.13) 

Proof.  Using  a  and  a  as  in  (2.2)  and  (2.12),  respectively,  we  have 

M  (0  =  max  max 

"P  ztQm  Kj<2p 

G/zi) 

from  (2.12) 

<  max  max 

ztQm  l«j«2p 

GfzS) 

from  (2.12) 

<  max  |e  (z;a)| 
*Qm 

=  E„(0 

implied  by  (2.7) 

from  (2.2) 

<  max  |e  (z;a)| 
*Qm 

from  (2.2) 

=  max  )  max 
ztQm  l  Kj<2p 

G^zia)  |  sec(tt/2p) 

implied  by  (2.7) 

=  Mnp(0  sec(n/2p)  . 

This  concludes  the  proof. 


Theorem  2.  Let  p>2  be  an  integer,  and  let  6S  =  n(j-l)/p,  j  =  1,2,  ....  2p.  Let 


^nD(f)  =  max  |e  (z;a)| 

where  the  complex  vector  atCn  is  any  vector  satisfying  (2.12).  Then 

(2.14) 

En(f)<^np(0<En(0sec(Tt/2p)  . 

(2.15) 

Proof.  Using  a  and  a  as  before,  we  have 

E  (f)  <  <f( f) 

n  np 

from  (2.2) 

=  max  |e  (z;&)| 

“Qm 

from  (2.14) 

<  max  |  max  G.(z;a)  ( sec(n/2p) 
z£Qm  (Kj<2p  ) 

implied  by  (2.7) 

<  max  max  G.(z;a)  sec(n/2p) 
ztQm  Kj<2p 

from  (2.12) 

<  max  |e  (z;a)|sec(n/2p) 

implied  by  (2.7) 

=  En(0  sec(n/2p). 
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This  concludes  the  proof. 


Corollary  2. 1.  Under  the  conditions  of  Theorem  2, 


M  (f)<  E  (f)<  6  (f)  . 

npv  n'  '  npv  ' 


(2.16) 


Proof.  Immediate. 


The  preceding  corollary  evidently  gives  excellent  upper  and  lower  bounds  on  the 
discrete  linear  Chebyshev  approximation  error  En(f),  and  these  bounds  are  readily 
available  after  the  numerical  computation  of  arCn  and  Mnp(f)  has  been  completed. 
We  point  out  that  the  above  two  theorems  substantially  generalize  results  in 
[3,  p.8541. 


Using  the  Maclaurin  series  for  sec  x  in  (2. 1 5)  gives  the  relative  discrepancy 

n  ^np(0  “  En(0  /  /-»  \  i  (  1  \ 

°< - — - «s<rc(*/2p>-l  +o(-).p-« o. 


(2.17) 


Note  that  this  upper  bound  on  the  relative  error  is  independent  of  f,  the  point  set 
Qm,  the  basis  functions  { hk } ,  and  n. 

We  will  now  explicitly  formulate  an  overdetermined  system  of  real  linear 
equations  to  be  solved  in  the  Chebyshev  norm  (to  be  defined)  which  is  equivalent  to 
solving  the  problem  (2.12).  Referring  to  the  choice  of  Oj’s  in  Lemma  2,  we  observe 
that  0p+j  =  7i  +  0j,  j  =  1,  .  .  .  ,  p,  and,  and  so  from  (2.11 ),  we  have 


Gp+j(z;a)  =  -Gj(z;a)  ,j  =  l . p 


Therefore,  we  may  rewrite  (2.12)  as 


M  (f)  =  min  max  |G.(z.;a)|  . 


a£Cn  l<i<m 


(2.18; 


Now,  breaking  the  following  quantities  into  their  real  and  imaginary  components 


we  may  write 


f(z)  =  u(z)  +  iv(z) 

hk(z)  =  rk(z)  +  isk(z)  ,  k  =  1,  .  .  .  ,  n, 
ak  =  bk  +  ick,  k  =  1,  .  .  .,  n, 


Rn(z;a)  =  u(z)  -  Z  b.r.(z)  +  Z  ck  s  (z) 

n  k=i  k  k  k=i 

n  n 

I„(z;a)  =  v(z)  -  2  bksk(z) "  ^  cv  rk<z)  • 


(2.19) 


(2.20) 


Using  (2. 1 1)  and  (2.20)  gives 
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G(7  :a)  -  mz^cosQ^  +  v(z()sin0j 

n 

-Z  bk  [rk(z|)  cos  +  sk(zt)  sin  a] 

n 

“Z  ck  [rk(zi)sin0j-sk(z|)cose].  (2.21) 

Note  that  Gj(zt;  a)  is  a  real  linear  equation  in  the  2n  variables  {bk}  and  {ck}, 
k  =  1,  .  .  .  ,  n,  and  that  all  the  coefficients  of  this  equation  are  computable  directly 
from  known  data. 

Define  the  mp  x  2n  real  matrix  B  in  the  partitioned  form 


B  = 


B.  D, 

B,  D, 


B  D 

m  m 


with  the  p  x  n  submatrices 

B,  =  [b<'>]  and  D,  =  [d«»J  ,  t  =  1 _ ,m 

whose  general  real  entries  are 


bjk  =  rk<z,>cos®j  +  W5"10, 

d\"=  rk(z()sinersk(zt)cose.] 
Also,  define  the  real  vector 

8  [s  1 1 »  •  *  •  '  *  *  ’ 

of  length  mp,  where 


j  =  1 . p;  k  =  1 . n. 

(2.22) 


’  gml’  '  '  ’  ’  gmpl 


(2.23) 


gtj  =  u(zt)  cos  6j  +  v(z()  sin  0.,  t  =  1 ,  .  .  .  ,  m;  j  =  1 . p. 

Finally,  define  the  real  vector 


x  =  [b,,  . 


•  ’  k„’  cr 


(2.24) 


of  length  2n.  With  this  notation  in  hand,  it  is  easily  seen  that  the  overdetermined 
system  of  mp  equations  in  2n  unknowns 


Bx  =  g 


(2.25) 
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has  a  residual  error  vector  defined  by 


g  -  Bx 

whose  mp  components  are  precisely  the  mp  real  numbers  Gj(2t;  a)  arranged  in  a 
special  order.  Therefore,  the  problem  (2.18)  can  be  solved  by  computing  a  solution 
to  the  overdetermined  linear  system  (2.25)  in  the  Chebyshev  norm;  i.e.,  the  largest 
magnitude  component  of  the  residual  vector  g-Bx  is  minimized  over  all  choices  of 
the  vector  x. 

This  equivalent  problem  in  linear  algebra  can,  in  principle,  be  solved  exactly  and 
in  a  finite  number  of  steps  using  linear  programming  methods  [1],  [3].  The  proof  of 
this  fact  is  the  content  of  the  following  self-contained  mathematical  result. 

Theorem  J.  Let  A  =  [a^]  be  a  real  r  x  s  matrix  with  r^s^l,  and  let  b  =  (b,,  .  .  .  ,  br) 
be  a  real  vector  of  length  r.  Let  a*  .  .  .  ,  a*+  ,,  co*  denote  a  solution  of  the  following 

primal  linear  program  in  the  s  +  2  real  variables  a{ . as+1,  w  w't*1  *'near 

constraints: 

Minimize:  co 

subject  to:  cr,  ^  0,  .  .  .  ,  a%  + ,  >  0,  co  >  0  , 


2i*kajk  +  «s,1A.  +  o>>bj,  j  =  l . r, 

S 

~£1Crkaik'"s+iAj  +  j  =  l . r’ 

where 

S 

A  =  -Z  a  j  =  l,  .  .  .  ,r  . 

J  k  =  i  JK 

Define 

Xk  =  °k~  as  +  t  ’  k  =  1.  •  •  •  .  s  . 

If 

S 

M  =  min  max  |b.-£  a  y  | 

Kj<r  J  k  =  1  JK 

with  the  minimum  taken  over  all  real  y ,,  .  .  .  ,  ys,  then 

M  =  co*  =  max  |b  a..x.| 
l<j<r  1  k=l 

where  x,,  .  .  .  ,  xs  are  given  by  (2.28). 


(2.26) 


(2.27) 

(2.28) 


(2.29) 


(2.30) 
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Proof.  We  first  prove  that  From  (2.28),  we  have  o*  = 

substituted  into  the  constraints  (2.26)  and  using  (2.27),  gives 


S 

I 


a.,  x. 

jk  k 


+  oj*  >  b.  ,  j  =  1,  .  . 


,  r  , 


and 


-I  a 


jk  Xk 


+  co*  >  -b.  ,j  =  1,  .  . 


,  r  . 


Clearly,  these  last  sets  of  inequalities  together  imply 


max 

l<j<r 


s 


(2.31) 


Hence,  from  (2.29),  M<oj*. 


We  next  prove  that  co*<M.  Let  x*,  .  .  .  ,  x*  denote  any  solution  of  (2.29).  Then 
we  may  write 


or,  without  absolute  values, 

br£,Vx*<M 

,  j  = . r  (2.32) 


Now,  define 


P  =  max  {0,-minx*}  (2.33) 

-  '  l«k<s  k 

=  Xk  +  *  i  ’  k  =  1 ,  .  .  .  ,  s  . 

Clearly,  the  s  +  2  real  numbers  /), . /3S4.,,  M  are  non-negative  by  construction. 

Furthermore,  substituting  x*  =  Pk  -  /Js  +  |,  k  =  1,  .  .  .  ,  s,  into  (2.32),  and  using 
(2.27),  gives 


b,  - 

Clearly  these  inequalities  show  that  the  numbers  ps+I,  M  satisfy  all  the 

constraints  (2.26).  Hence,  it  must  be  that  oj*<M.  Since  we  have  already  established 
that  M<oj*,  we  conclude  that  M  =  cu*.  Hence  the  inequality  (2.31)  must  actually  be 
an  equality  in  light  of  the  definition  (2.29).  This  completes  the  proof. 
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Theorem  3  does  not  require  that  the  solutions  of  either  the  linear  program  or  of 
(2.29)  be  unique.  Theorem  3  states  only  that  from  a  given  solution  of  the  linear 
program,  we  may  construct  a  solution  of  (2.29)  using  (2.28).  Conversely,  it  is  easy 
to  see  that  any  solution  of  (2.29)  can  be  used  to  construct  a  solution  of  the  linear 
program  using  (2.33). 

An  excellent  algorithm,  which  we  will  refer  to  as  ACM  495,  is  available  in  the 
literature  [1]  for  solving  the  linear  program  of  Theorem  3.  It  requires  as  input  only 
the  overdetermined  system  of  equations  Ax  =  b.  The  linear  program  is  then  set  up 
and  solved  by  the  algorithm,  so  that  knowledge  of  linear  programming  techniques  is 
not  necessary  to  use  the  algorithm  in  practice.  The  computational  procedure,  in¬ 
ternal  to  the  algorithm,  actually  solves  the  dual  of  the  above  primal  linear  program 
using  a  modification  of  the  simplex  method.  The  dual  formulation  of  this  problem 
is  available  in  [5,  p.  296],  We  will  not  discuss  the  details  of  the  linear  programming 
technique  further  in  this  report. 

A  very  simple  modification  [3,  p.  863]  of  ACM  495  yields  an  algorithm  for 
solving  any  real  overdetermined  system  of  linear  equations  in  the  Chebyshev  norm 
subject  to  the  additional  constraints  that  all  the  residuals  be  non-negative.  For  A 
and  b  as  in  Theorem  3,  this  problem  takes  the  form 

minimize  max  (b  a...  x.  ^  (2.34) 

x, . x5  l<3<rv  J  k  =  l  Jk 

subject  to  the  r  constraints 

b:-Z  aikx.  >0,  J  =  1 . r-  (2-35) 

The  solution  x,,  .  .  .  ,  xs  returned  by  this  modified  algorithm  is  correct,  even  though 
the  residuals  returned  may  be  in  error.  The  correct  residuals,  if  desired,  must  be 
calculated  directly  from  the  solution.  Alternatively,  if  the  residuals  are  required  to 
be  non-positive, .then  the  same  modified  algorithm  will  work  with  A  and  b  replaced 
by  -A  and  -b,  respectively. 

Requiring  non-negative  residuals  in  the  overdetermined  system  (2.25)  has  in¬ 
teresting  geometrical  interpretations.  For  example,  if  we  take  p  =  2  in  Lemma  2, 
then  0,  =  Oand02  =  rc/2.  Thus,  from  (2.1 1),  G,(z;a)  and  G,(z;a)  are  merely  the  real 
and  imaginary  parts  of  the  complex  error  en(z;a).  Thus,  the  2m  components  of  the 
residual  vector  g-Bx  are  precisely  the  real  and  imaginary  parts  of  en(z;a)  evaluated 
at  all  m  data  points.  Therefore,  if  the  system  (2.25)  is  required  to  have  non-negative 
residuals,  we  have  forced  the  error  curve  to  lie  entirely  in  the  first  quadrant  of  the 
complex  plane.  More  generally,  we  may  always  constrain  en(z;a)  to  lie  in  a  given 
convex  wedge  shaped  sector  W  of  the  complex  plane  with  vertex  at  the  origin,  by 
making  different,  but  appropriate,  choices  of  the  angles  0!  and  0,.  Further  ex¬ 
ploration  of  this  idea  shows  that  upper  and  lower  bounds  for  the  error  Wn(f), 
defined  by 
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W  (f)  s  min  max  |e  (z;a)| 
act"  «Qm  ' 

subject  to:  en(z;a)  z# ,  ztQm 

where  en(z;a)  is  given  by  (2.1),  can  be  obtained  in  terms  of  the  error  Wnf>( f),  defined 
by 


W  (0  s  min  max  max  G.(z;a) 

"P  a tCn  z£Qm  j  =  i p 

subject  to:  G^fzja)  >  0,  Z£Qm,j  =  l . p, 

where  G^zta)  is  given  by  (2.21),  but  for  a  different  set  of  angles  0(,  .  .  .  ,  0  .  The 
quantity  W  (f)  may  be  computed  by  solving  the  linear  system  analogous  to  (2.25) 
in  the  Chebyshev  norm  with  the  constraint  of  non-negative  residuals.  This  approach 
is  especially  effective  when  the  vertex  angle  of  the  wedge  does  not  exceed  n/2.  This 
topic  is  hot  pursued  further  in  this  report. 

Suppose,  finally,  that  the  complex  solution  vector  aeCn  of  problem  (2.12)  is 
required  to  be  strictly  real,  while  f  and  {hk}  are  complex.  Then,  in  the  vector  x  of 
(2.24),  c,  =  .  .  .  =  cn  =  0.  Thus,  the  overdetermined  system  Bx  =  g  of  mp 
equations  in  2n  unknowns  can  be  replaced  by  a  smaller  system  &x  =  g  of  mp 
equations  in  only  n  unknowns,  where  the  mp  x  n  real  matrix  £  is  defined  in 
partitioned  form  by 


B, 


B 


(2.36) 


where  the  pxn  submatrices  B,,  .  .  .  ,  B  are  unchanged  from  (2.22),  and  the  real 
vector  x  =  [b,,  .  .  .  ,  bn]T.  A  solution  of  §£  =  g  in  the  Chebyshev  norm  can  be 
computed  using  linear  programming  and  algorithm  ACM  495  as  before. 


III.  Numerical  Examples  and  Efficiency  of  Approach 


We  illustrate  the  procedure  of  the  preceding  section  by  approximating  the 
complex  function  f(x)  =  exp  (i3x)  by  a  weighted  sum  of  the  basis  functions  1, 
exp(ix),  exp(i2x).  That  is,  we  seek  to  minimize  the  magnitude  of  the  complex  error 
curve 

e,(x)  =  exp(i3x)  -  2  a.  exp(i(k-l)x)  (3.1) 

k  =  i 

over  interval  [0,  n/4],  by  choice  of  a,,  a,,  a3,  by  solving  the  problem  M  (0  of 
(2.12).  Two  cases  are  of  interest;  in  the  first,  the  coefficients  {ak}j  are  restricted  to 


II 


be  real,  whereas  in  the  second,  these  coefficients  can  be  complex.  1  he  number  m.  ol 
equtspaced  x-xalues  at  which  (  V I )  is  sampled,  is  taken  to  be  either  II,  101,  or  1001 , 
thereby  ensuring  that  the  smaller  sample  sizes  are  subsets  ol  the  larger  si/es.  lire 
value  ol' p.  which  is  hall  the  number  ol  phase-shifted  values  ol  (.V I )  employed  in  lire 
error  minimization.  is  taken  to  be  2,  6,  IS,  54,  again  ensuring  the  subset  behav tor  of 
the  smaller  si/e  cases.  Note  that  p  and  the  phase  shifts,  9  | ,  are  as  given  in  Theorem 
1  m  Section  II 

The  optimum  real  coefficients  in  (3.1)  for  the  problem  Mnp(f)  are  given  in  table  1 
for  these  choices  of  m  and  p,  and  a  plot  of  the  magnitude  of  the  error  for  several 
representative  cases  is  given  in  figure  I .  The  best  approximation  of  all  cases  con¬ 
sidered  is  afforded  by  m  =  1001 .  p  =  54*  and  its  error  curve  is  plotted  as  a  solid  line; 
its  maximum  error  is  .1078,  which  is  realized  at  two  points  in  the  interval  |0,  n/4). 
The  cases  for  smaller  m  (less  sampling  of  the  abscissa)  and  smaller  p  (less  sampling 
of  the  phase  of  the  complex  error)  are  poorer;  for  example,  the  maximum  error  for 
m  =  1 1,  p  =  2  is  .1 184,  realized  at  only  one  point,  namely  x  =  n  4. 


Table  1.  Coefficients  for  the  Real  Weight  Case 


m 

P 

at 

a2 

a? 

11 

2 

.936738 

-2.443144 

2.518388 

6 

.828404 

-2.280319 

2.396455 

18 

.858547 

-2.321885 

2.425096 

54 

.844146 

-2.301461 

2.410611 

101 

2 

.936781 

-2.443223 

2.518458 

6 

.831314 

-2.284548 

2.399525 

18 

.865131 

-2.331446 

2.432033 

54 

.853823 

-2.315301 

2.420506 

1001 

2 

.936785 

-2.443232 

2.518466 

6 

.831237 

-2.284448 

2.399461 

18 

.865213 

-2.331571 

2.432127 

54 

.853443 

-2.314772 

2.420138 

We  have  not  plotted  the  other  error  curves  with  real  coefficients  tor  m  =  101  and 
1001 ,  because  they  are  indistinguishable  from  figure  1 ,  as  a  perusal  of  table  I  shows. 
For  example,  the  coefficients  for  m  =  1 1 .  p  =  2  are  very  close  to  those  for  m  =  101 . 
p  =  2  and  m  =  1001 ,  p  =  2.  Thus,  our  sampling  in  x  is  already  “fine  enough”  at 
m  =  1 1 .  However,  there  is  a  significant  change  in  the  coefficients  as  p  is  varied,  for  a 
fixed  value  of  m;  that  is,  p  =  2  yields  very  coarse  phase-sampling  of  the  error  curve 
and  should  definitely  be  made  larger. 

The  Chebyshev  error  curve  (m  =  1001,  p  =  54)  in  figure  1  realizes  its  maximum 
value  at  only  n-1  points,  rather  than  at  n  +  1  points,  where  n  =  3  is  the  number  of 
coefficients  for  this  example.  This  is  probably  related  to  the  fact  that  we  have 
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Figure  I.  Frror  Curves  fur  Real  C  oefficients;  m  =  1 1 

minimized  both  the  real  and  imaginary  parts  of  the  complex  error,  but  have  allowed 
ourselves  to  use  only  real  coefficients. 

The  solution  of  the  problem  M|>p(f)  for  complex  weights  is  given  in  table  2  for  the 
same  choices  of  m  and  p  as  above.  Again,  the  change  in  coefficient  values  is  more 
marked  w tth  p  than  with  m.  Magnitude-error  curves  for  m  =  1 1  and  101  are  given  in 
figuies  2  and  3,  respectively;  the  curves  for  m  =  1001  are  indistinguishable  from 
those  for  m  =  101  and  are  not  presented. 


Table  2.  Coefficients  for  the  Complex  Weight  Case 


m 

P 

Re(a,) 

lm(a, ) 

Re(a:) 

Im(a:) 

Re(a,) 

Im(a,) 

11 

2 

.364737 

.954343 

-2.021670 

-2.119639 

2.669023 

1.153207 

6 

.378045 

.907888 

-2.016657 

-2.018598 

2.648834 

1.100488 

18 

.373079 

.898715 

-2.003032 

-2.003205 

2.639992 

1.094451 

54 

.371586 

.896504 

-1.999352 

-1.999473 

2.637788 

1 .092947 

101 

2 

.362962 

.953469 

-2.018255 

-2.119960 

2.667544 

1.154238 

6 

.376532 

.904026 

-2.012095 

-2.014055 

2.646131 

1.099461 

18 

.370549 

.893500 

-1.995913 

-1.997062 

2.635782 

1.093144 

54 

.368950 

.890017 

-1.991172 

-1.99)196 

2.632622 

1 .090777 

1001 

2 

.362947 

.953499 

-2.018253 

-2.120028 

2.667560 

1.154275 

6 

.376502 

.903926 

-2.011979 

-2.0139)4 

2.646047 

1.099417 

18 

.370711 

.893848 

-1.996440 

-1.997545 

2.636145 

1.093278 

54 

.369179 

.890566 

-1.991954 

-1.991974 

2.633175 

1.091006 
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The  Chebyshev  error  eurve  (m  =  1001,  p  =  54)  is  now  symmetric  about  the 
midpoint  of  the  interval  of  interest  and  has  four  equal  error  peaks  of  value  .0147. 
This  error  is  7.3  times  smaller  than  that  for  the  real  coefficient  case.  Also  the 
number  of  equal  error  peaks  now  equals  I  plus  the  number  of  coefficients;  whether 
this  property  holds  generally  is  not  known. 

Upper  and  lower  bounds  on  the  discrete  Chebyshev  error  En(f)  for  the  real  and 
complex  coefficient  cases  are  given  in  table  3.  These  bounds  are  precisely  those 
presented  in  (2.16).  They  correspond  to  sampling  the  complex  error  (3.1 )  both  in  the 
abscissa  x  and  in  the  phase  of  e,(x).  The  lower  bounds  monotonically  increase  with 
increasing  m  or  p.  The  upper  bounds  decrease  with  increasing  p,  but  increase  with 
increasing  m.  All  these  trends  follow  from  the  fact  that  smaller  sample  sizes  are 
subsets  of  the  larger  sizes. 

Table  3.  Bounds  on  the  Discrete  Chebyshev  Error  En(f) 


m 

P 

Real  Coefficients 

Lower  Bound  Upper  Bound 

Complex  Coefficients 
Lower  Bound  Upper  Bound 

11 

2 

.083718 

.118396 

.012089 

.017097 

6 

.105074 

.108780 

.013963 

.014456 

18 

.107307 

.107717 

.014143 

.014197 

54 

.107612 

.107658 

.014168 

.014174 

101 

2 

.083731 

.118414 

.012252 

.017328 

6 

.105192 

.108893 

.014436 

.014946 

18 

.107556 

.107967 

.014677 

.014733 

54 

.107767 

.107813 

.014703 

.014709 

1001 

2 

.083734 

.113418 

.012255 

.017331 

6 

.105191 

.108901 

.014440 

.014950 

18 

.107565 

.107976 

.014679 

.014735 

54 

.107775 

.107821 

.014704 

.014712 

However,  the  maximum  magnitude  error,  evaluated  over  the  continuum  of  x- 
values  in  the  interval  [0,  n/4]  (actually  computed  on  a  dense  discrete  sampling 
space),  obeys  none  of  these  monotonic  relations,  as  table  4  demonstrates.  For 
example,  the  maximum  error  in  the  real  case  for  m  -  1 1,  p  =  18  is  less  than  that  for 
m  =  ll,  p  =  54.  Also,  the  maximum  error  in  the  complex  case  for  m  =  ll.  p  =  6  is 
greater  than  that  for  m  =  101,  p  =  6.  The  reason  for  this  behavior  is  that  we  have 
minimized  a  discrete  approximation  to  our  problem  of  interest,  sampling  both  in  the 
abscissa  x  and  in  the  phase  values  of  the  complex  error.  However,  the  numerical 
discrepancies  are  small,  as  they  must  be  for  reasonably  fine  sampling  in  both 
variables.  (A  recursive  gradient  procedure  could  be  used  with  any  of  these  coef¬ 
ficient  sets  to  improve  the  final  maximum  magnitude-error  if  desired.) 

The  FORTRAN  program  listing  in  the  Appendix  is  the  exact  code  used  to 
generate  the  complex  weights  in  example  (3.1)  for  m  =  101  and  p  =  6.  The  imbedded 
comments  should  enable  anyone  seeking  to  use  and  understand  the  code  to  do  so. 
Further  remarks  are  given  in  the  discussion  in  Section  VI. 
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Table  4.  Maximum  Magnitude  Error, 
Computed  Over  2001  Equispaced  Points  in  |0,  rt/4| 


m 

P 

Real  Coefficients 

Complex  Coefficients 

11 

2 

.118396 

.017097 

6 

.108780 

.015142 

18 

.107890 

.015004 

54 

.107983 

.015005 

101 

2 

.118415 

.017329 

6 

.108893 

.014946 

18 

.107967 

.014733 

54 

.107813 

.014711 

1001 

2 

.118417 

.017331 

6 

.108902 

.014950 

18 

.107976 

.014735 

54 

.107821 

.014712 

Ef fit’s-  and  timing  estimates  for  actual  calculation  of  complex  Chebyshev 

a|'orssitinu'1jns  by  the  method  of  this  report  is  an  important  consideration  in  some 
applications.  If  we  define  an  operation  as  consisting  of  a  multiplication  followed  by 
an  addition,  then  it  is  known  [6]  that  the  number  of  operations  per  simplex  iteration 
requi's  J  by  algorithm  ACM  495  [1]  is  exactly  the  number  of  equations  times  the 
number  of  unknowns.  In  our  case,  the  number  of  equations  is  mp,  and  the  number 
of  unknowns  is  2n  if  the  coefficients  are  complex,  or  n  if  the  coefficients  are 
required  to  be  real.  Thus,  the  operation  count  per  iteration  is  either  2nmp  or  nmp. 
The  number  of  iterations  required  is  difficult  to  estimate,  since  it  depends  on  the 
particular  problem.  However,  in  randomly  generated  problems,  it  has  been  ob¬ 
served  [6]  that  the  number  of  iterations,  I,  is  approximately  the  number  of 
unknowns  times  some  small  constant  c,  where  usually  l<c<3.  (Similar  estimates 
have  been  observed  [7,p.  160],  [8]  in  more  general  linear  programs  as  well.)  Thus,  in 
our  case,  I  =  2cn  if  the  coefficients  are  complex  and  I  =  cn  if  they  are  real. 

The  CPU  time  should  be  proportional  to  the  total  operation  count,  which  equals 
the  product  of  the  number  of  iterations  and  the  number  of  operations  per  iteration. 
That  is,  we  expect  the  CPU  time  to  be  proportional  to  n:mp.  For  the  particular 
example  here,  however,  we  obtain  an  excellent  fit  to  the  limited  data  in  table  5  with 

CPU  time  (msec)  =  .128  n1  13  m1  18  p1  ,8> 

where  n  =  6  if  the  coefficients  are  complex,  and  n  =  3  if  they  are  real.  This  fit  was 
obtained  by  letting  the  exponents  of  n,  m,  and  p  vary  separately.  Other  examples, 
however,  lead  us  to  anticipate  that,  more  generally, 

CPU  time  oc  n;(mp)i  ;, 

with  a  proportionality  factor  of  the  order  of  .01-  03  msec,  w  here  n  is  either  twice  the 
number  of  approximation  coefficients  if  the  coefficients  are  complex,  or  exactly  the 
number  of  coefficients  if  they  are  required  to  be  real. 
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Table  5.  Number  of  Simplex  Iterations  and  CPU  Time 


m 

P 

Real  Coefficients 
Simplex  CPU  (s) 

Complex  Coefficients 
Simplex  CPU  (s) 

11 

2 

6 

.02 

10 

.05 

6 

8 

.08 

15 

.16 

18 

11 

.23 

21 

.58 

54 

13 

.81 

27 

2.25 

101 

2 

7 

.25 

10 

.40 

6 

9 

.73 

17 

1.60 

18 

13 

2.65 

21 

5.78 

54 

15 

11.39 

28 

24.27 

1001 

2 

9 

3.05 

13 

5.00 

6 

10 

10.34 

17 

19.38 

18 

13 

48.16 

24 

105.47 

54 

16 

170.52 

28 

359.20 

The  CPU  time  estimates  apply,  of  course,  only  to  the  DEC  VAX  1 1/780  com¬ 
puter  on  which  the  calculations  were  performed.  The  virtual  memory  feature  of  this 
system  allows  very  large  problems  to  be  solved;  however,  for  sufficiently  large 
problems,  the  system  overhead  incurred  (page  faulting,  and  so  on)  may  significantly 
and  adversely  affect  these  estimates. 

IV.  Application  to  Array  Design  with  a  Constraint 

Consider  a  linear  antenna  array  with  N  elements  located  at  arbitrary  fixed 
positions  {xk}^  receiving  a  plane  wave  arrival  of  wavelength  A  from  direction  9a, 
--y  <0a<-y>  relative  to  a  normal  to  the  array.  If  the  array  is  steered  to  look  in 
direction  6,,  < 6 then  the  complex  transfer  function  of  the  beamformer  is 

given  by 

N 

T(u)  =  2  wkexp(-idku)  , 

k  =  l  (4.1) 

where  { wk }*are  the  element  weights,  and 

dk=2nxk/A  for  l<k<N  , 

u  =  sin  0a  -  sin  0j  . 

Observe  that  the  total  range  of  u  depends  on  the  look  direction  0| ;  for  example,  if 
=  0,  then  the  range  of  u  is  the  closed  interval  [-1,1],  The  peak  response  of  T(u) 
should  occur  at  u  =0,  so  we  normalize  (without  loss  of  generality)  according  to 

N 

T(0)  =1=2  W 

k  =  i 
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To  realize  small  sidelobes,  we  must  minimize  !T(u)|  for  all  u  values  in  some  subset 
U  of  the  total  range  of  u.  f  or  example,  if  9(  =  0,  the  total  range  of  u  is  \- 1.1],  and 
U  could  be  the  union  of  intervals  [-1,-uJ  and  [u„,l],  where  u()>0  is  chosen  small 
relative  to  1.  For  the  special  case  of  real  weights  {wk},  since  from  (4.1), 
T(-u)  =  T*(u),  we  can  confine  attention  to  U  =  [u„  ,  1].  The  normalization  con¬ 
straint  is  most  easily  accounted  for  by  solving  for  w  N  and  eliminating  it;  we  obtain 
then 

N-l 

T(u)  =  exp(-idNu)-2  wk  (exp(-idNu)  -  exp(-idku))  .  (4.2) 

k  =  I 

This  problem  now  fits  the  framework  of  (2.1)  if  we  identify 

z  =  u 
n  =  N-l 
en(z)  =  T(u) 
f(z)  =  exp(-idNu) 

ak  =  wk 

hk(z)  =  exp(-idNu)  -  exp(-idku) 

Qm  =  finite  subset  of  U  .  (4.3) 

There  has  been  no  statement,  thus  far,  as  to  the  real  or  complex  nature  of  the 
weights  {wk}.  This  distinction  depends  upon  the  application  and  the  capability  of 
the  beamformer.  Both  cases  fit  the  above  framework;  the  only  difference  is  that  the 
number  of  unknowns  to  be  solved  for  will  be  twice  as  large  for  the  complex  weights 
as  for  the  real  weights. 

If  the  array  is  half-wavelength  equispaced,  then  the  computed  element  weights 
will  be  identical  to  the  classical  Dolph-Chebyshev  weights  and  can,  in  this  instance, 
be  computed  analytically.  The  general  case  of  arbitrary  spacings,  however,  cannot 
be  computed  analytically,  yet  the  algorithm  presented  in  this  report  can  always  be 
applied. 

In  the  remainder  of  this  section,  we  presume  that  the  elements  are  equispaced  at 
half-wavelength.  Then  xk  =  kA/2  and  (4.1)  becomes 

N 

T(u)  =  2  wt  exp(-inku)  .  (4.4) 

k  =  I 

Observe  now  that  T(u)  in  (4.4)  has  period  2  in  u,  regardless  of  whether  the  weights 
{wfc}  are  real  or  complex,  or  whether  some  elements  have  failed,  i.e.,  zero  weight 
values.  This  means  that  we  can  study  and  control  T(u)  in  (4.4)  over  any  convenient 
u-interval  of  length  2,  and  need  not  confine  our  investigation  to  [-1,11.  In  par¬ 
ticular,  we  concentrate  on  the  u-interval  [0,2]  in  the  following. 

As  an  illustration  of  the  capability  of  the  minimization  technique  of  this  report,  a 
50-element  half-wavelength  equispaced  linear  array  was  initially  designed  for  peak 
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sidelobes  of  -30  dB  relative  to  the  main  peak.  This  is  of  course  a  standard  Dolpli- 
Chebyshev  case,  and  gives  -30  dB  sidelobes  throughout  the  u-range  (u(l,  2-u„], 
where  u„  =  .0538117.*  Then  10°’o  of  the  elements  were  randomly  eliminated  from 
the  array,  but  the  remaining  weights  were  unchanged;  this  corresponds  to  5  elements 
failing  in  the  array.  The  relative  response  of  this  particular  array,  with  elements  7, 
22,  40,  43,  50  failed,  is  illustrated  in  figure  4.  The  peak  sidelobe  has  increased  from 
-30  dB  to  -21.58  dB,  a  degradation  of  8.4  dB,  and  there  is  a  large  variety  of  dif¬ 
ferent  size  peaks. 


0  .25  .5  .75  1  1.25  1.5  1.75 

u  =  sin  9g  —  sinej 


Figure  4.  Relative  Pattern  for  5  Elements  Failed 

When  our  method  with  p  =  2  and  m=  251  equispaced  points  in  [u0,  2-u0]  is  ap¬ 
plied  to  this  defective  array,  and  the  remaining  45  elements  are  weighted  with  real 
coefficients,  subject  to  the  constraints  that  the  mainlobe  width  be  the  same  as  the 
ideal  50-element  array,  and  that  the  steering  range  in  u  be  the  same,  the  resultant 
array  pattern  is  displayed  in  figure  5.  The  peak  sidelobe  is  now  -23.62  dB,  an  im¬ 
provement  of  2.04  dB  over  figure  4;  however,  there  is  still  a  significant  variation  in 
the  values  of  the  sidelobes,  due  to  an  insufficient  number  of  phase  controls,  namely 
only  p  =  2. 

When  we  increase  the  parameter  values  to  p  =  8,  m  =  501,  the  resultant  best  real 
weights  are  displayed  graphically  in  figure  6  and  the  corresponding  array  pattern  is 
given  in  figure  7.  The  gaps  in  figure  6  at  locations  7,  22,  40,  43,  50  correspond  to 
zero  weighting  at  the  failed  elements.  The  general  character  of  the  weights  is  a  bell- 


*For  an  N-element  array  and  -t  dB  peak  sidelobes.  we  have  u0=  :2/n)  arc  cos(l/zQ)  where 
2z0  =  |r  +■  ] 1  M  +|r-»Jr^T  ]1  M,  r=l0*  :o,andM  =  N-l. 
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0  .25  .5  .75  1  1.25  1.5  1.75  2 

u  =  sin  Sa  -  sin  0j 

Figure  5.  Relative  Pattern  for  p  =  2,  m  =  251,  Real  Weights 


Figure  6.  Best  Real  Weights  for  p  =  8.  m  =  501 
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0  .25  .5  75  1  1.25  1.5  1.75  2 


u  =  sin8a  -  sin  8j 

Figure  7.  Relative  Pattern  for  p  =  8.  m  =  501,  Real  Weights 

shaped  one  of  all  positive  numbers,  but  there  is  significant  fluctuation  in  the  actual 
weight  values,  of  the  order  of  lO^o.  The  pattern  in  figure  7  has  a  peak  sidelobe  of 
-25.20  dB,  an  improvement  of  3.62  dB  over  figure  4,  but  still  4.80  dB  poorer  than 
the  ideal  50-element  array. 


When  the  weights  were  allowed  to  be  complex,  and  the  maximum  sidelobe 
minimized  in  the  same  steering  range  [u0,  2-u0]  for  p  =  2  and  m  =  501  equispaced 
points  in  [u0,  2-u0],  the  best  complex  weights  turned  out  to  be  virtually  pure  real, 
and  the  corresponding  pattern  was  almost  identical  to  figure  5.  A  much  improved 
pattern  for  complex  weights  was  achieved  when  we  took  p  =  8,  m  =  501 ;  in  fact,  the 
best  complex  weights  were  real  (within  lO-6  relative  error)  and  the  pattern  was  the 
same  as  figure  7.  Although  we  had  anticipated  a  better  pattern  for  the  complex 
weight  case  than  for  the  real  weights,  that  did  not  materialize;  the  best  complex 
weights  for  this  equispaced  linear  array  with  5  missing  elements  were  real.  The 
reason  for  this  behavior  is  unknown,  but  it  is  an  encouraging  result  from  the  array 
design  viewpoint,  for  it  indicates  that  there  is  no  need  to  allow  phasing  at  the  in¬ 
dividual  elements;  gain  alone  will  achieve  all  the  sidelobe  reduction  that  can  be 
achieved.  This  conclusion  is  drawn  only  for  the  half-wavelength  equi-spaced  line 
array  with  omnidirectional  element  response. 

The  use  of  linear  programming  to  design  antenna  arrays  is  not  entirely  new.  In  [9] 
and  [10],  linear  programming  was  used  to  synthesize  desired  complex  transfer 
functions  to  within  3  dB  of  the  best  possible  sidelobe  level.  Their  method 
corresponds  to  taking  p  =  2  in  the  method  presented  in  this  report. 
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Tin.'  compulation  of  the  real  weight'  ol  figure  6  (w here  p  2.  m  -  25 1 .  and  n  =  44) 
and  of  figure  7  (where  p  -  8,  in  501,  and  n  44)  required  1.2  minutes  205 
iterations  and  38.4  minutes  402  iterations,  re-pe i v els .  On  the  other  hand,  when 
the  weights  were  allowed  to  be  complex  (replacing  n  =44  bv  n  =88,  but  leaving  p 
and  m  unchanged  in  both  cases),  the  computations  required  7.0  minutes  657 
iterations  and  179  minutes  '  262  iterations,  respectively  .  The  two  ol  these  tour  cases 
requiring  the  smallest  CPU  dines  encountered  almost  no  s\ stein  overhead  due  to 
program  size.  However,  the  two  cases  requiring  the  largest  CPU  times  encountered 
very  significant  system  overhead  because  their  large  memory  requirements  caused 
significant  usage  of  the  virtual  memory  feature  of  the  DEC  VAX  1 1/780.  The  38.4 
minute  case  required  over  3.6  million  page  faults,  while  the  179  minute  case 
required  over  11  million  page  faults.  It  is  important  to  bear  in  mind  that  the  DEC 
VAX  1 1/780  is  essentially  a  mini-computer,  and  that  without  virtual  memory,  only 
the  largest  mainframe  computers  could  have  solved  either  of  these  two  problems. 


V.  Efforts  to  Extend  the  Method 

Our  basic  problem  is  to  minimize  the  maximum  magnitude  of  complex  error 

en(z)  =  f(z)  -  Z  (  ak  hk(z)  (5.1) 

over  a  continuum  of  values  of  z,  when  f,  { hk } ,  and  {ak}  are  complex.  We  im¬ 
mediately  approximate  this  desired  problem  by  discretizing  the  z  variable  to  a  finite 
number  of  values,  in  order  to  make  the  problem  computable.  Furthermore,  at  any  z 
value  of  interest,  we  additionally  discretize  the  number  of  phase  errors  we  are 
willing  to  consider.  To  be  specific,  since  the  algorithm  in  [1]  applies  only  to  real 
quantities,  we  consider  the  “projection”  of  a  rotated  version  of  the  complex  error: 

PU.H')  =  Re{exp(ivF)  en(z):  ,  (5.2) 

Then,  since  the  argument  of  complex  error  (5.1)  is  unknown  a  priori ,  we  let  V  take 
on  a  finite  set  of  values  spaced  over  any  rt  radian  interval,  and  minimize  the 
magnitude  of  projection  (5.2)  over  all  these  selected  V  values.  This  is  equivalent  to 
the  method  of  Section  II. 

In  an  effort  to  eliminate  this  second  discretization  process  in  V,  a  perturbation 
method  was  put  forth  in  [3]  that  claimed  guaranteed  convergence  to  the  optimum 
weights,  for  any  given  finite  discrete  set  of  z-values.  When  applied  to  the  examples 
in  [31,  the  proposed  perturbation  technique  did  indeed  converge.  However,  when 
applied  to  the  following  example,  of  approximation  of  exp(i3x)  by  the  three  basis 
functions  1,  exp(ix),  exp(i2x),  over  100  equispaced  points  in  the  domain  [0,  n/4]  in 
x,  it  sometimes  failed  to  converge,  depending  on  the  initial  weights  employed.  The 
reason  for  this  failure  is  that  the  “direction  of  the  minimum”  furnished  by  the 
perturbation  is  often  totally  irrelevant,  and  the  best  scale  factor  to  apply  to  this 
perturbation  is  very  small.  Thus  there  occurs  a  small  random  meander  in  the 
coefficient  space,  and  occasional  convergence  to  a  non-optimum  point. 
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A  modification  of  this  technique  was  attempted  wherein  the  magnitude  of  the 
perturbation  was  bounded.  Although  this  improved  the  situation  somewhat, 
convergence  to  the  optimum  was  not  always  obtained. 

It  was  thought  that  this  meander  in  coefficient  space  might  be  eliminated  by 
tracking  the  exact  z-values  at  which  (5.1)  is  a  maximum.  Recall  that  in  the  real  case 
discussed  in  the  Introduction,  convergence  to  the  absolute  optimum  over  a  con¬ 
tinuum  of  real  z-values  was  achieved  in  a  practical  example  by  re-evaluating  the  z- 
points  of  maximum  error  and  using  these  in  a  recursive  approach.  When  this  idea 
was  extended  to  the  two  continuous  variables  z,  V  in  (5.2),  and  only  the  2n  +  1 
largest  error  points  were  retained,  convergence  was  not  obtained.  When  however, 
the  single  “point”  of  a  maximum,  i.e.,  a  pair  of  values  (zk,  Vk),  was  replaced  by  a 
“patch”,  i.e.,  a  set  of  values  {(zkp,  4,kp)}  covering  the  maximum  point  (zk,  4>k),  the 
convergence  to  the  absolute  optimum  for  the  examples  considered  was  apparently 
achieved.  The  patch  width  in  was  of  the  order  of  a  degree  in  most  cases.  The 
problem  with  this  latter  modification  is  that  a  large  number  of  computations  of  the 
error  function  and  its  derivative  must  be  evaluated,  and  the  improvement  over  the 
method  of  Section  II  is  insignificant  when  p  there  is  large. 

If  the  final  error  in  (5.2),  after  application  of  the  method  of  Section  II,  is 
inadequate,  due  to  inadequate  sampling  in  z  and/or  V,  it  is  possible,  for  a  given 
coefficient  set  {ak},  to  locate  the  point  (zm,  Vm)  at  which  (5.2)  is  largest,  and  then 
use  a  gradient  approach  to  decrease  this  maximum  error  at  (zm,  VJ.  Of  course,  the 
particular  point  of  maximum  will  jump  around  as  the  set  {ak}  is  perturbed; 
nevertheless,  the  technique  does  converge  (although  slowly)  and  does  lead  t  >  smaller 
errors  at  the  maximum  of  (5.2)  in  a  continuum  for  z  and  ML 


VI.  Discussion  and  Summary 

It  has  been  observed  that  two  of  the  locations  of  maximum  magnitude  error  often 
occur  at  the  endpoints,  if  the  specified  domain  in  (1.2)  is  a  real  interval;  for 
example,  see  figures  2  and  3.  (The  example  of  real  coefficients  in  figure  1  had  one  of 
the  maximum  error  points  at  an  endpoint,  but  not  the  other.  However,  if  we  had 
specified  domain  [-n/4,  ti/4]  in  that  example,  we  would  have  observed  four  peak- 
error  points,  two  of  which  would  have  been  at  endpoints,  due  to  the  conjugate 
property  of  the  desired  function  and  the  basis  functions.)  Since  the  endpoints  may 
be  the  only  ones  we  can  anticipate  a  priori  and  specify  as  locations  of  maximum 
error,  an  obviously  useful  procedure  is  to  use  more  values  of  phase  shift  V  in  (5.2) 
(alternatively,  the  angles  {0  }  in  Lemma  2)  at  the  endpoints  than  in  the  interior,  so 
as  to  better  control  these  very-likely  locations  of  maximum  error.  For  example,  we 
might  use  p  =  6  in  the  interior  of  a  specified  real  interval  domain  of  z  and  use  p  =  12 
or  20  at  the  two  endpoints.  This  does  not  add  greatly  to  the  total  computation,  since 
there  are  generally  far  more  interior  points  than  (two)  endpoints.  The  program  in 
the  Appendix  may  be  readily  used  with  different  values  of  p  at  different  data  points 
by  exploiting  the  INDEX  array  in  the  user-supplied  subroutine  named  ZPHASE. 


The  p  different  phase  shifts  V  selected  in  (5.2)  have  been  chosen  here  to  be  equally- 
spaced  over  a  180°  span  (along  with  their  180°  mates).  This  is  the  most  reasonable 
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selection  in  the  absence  of  a  priori  knowledge  of  the  complex  error,  its  magnitude, 
and  phase  because  it  gives  the  best  upper  bound  (2.7)  in  Lemma  2  of  any  set  of 
phases.  However,  one  could  select  any  value  of  4*  to  investigate  the  error;  for 
example,  different  sets  of  values  of  4*  could  be  used  at  various  values  of  abscissa  z. 
The  program  in  the  Appendix  may  be  used  with  any  desired  set  of  phases  at  any,  or 
all,  of  the  data  points  simply  by  altering  the  user-supplied  subroutine  named 
ZPHASE. 

The  potential  for  significant  round-off  error  accumulation  is  always  present  in 
linear  Chebyshev  complex  function  approximation.  For  example,  in  approximating 
f(x)  =  cos  (12x)  +  i  sin  (3x)  by  a  complex  linear  combination  of  the  12  basis 
functions  1,  exp(ix),  ....  exp  (il  lx)  on  the  interval  [0,  rr/4],  the  complex  coef¬ 
ficients  of  best  approximation  were  observed  to  be  large  in  magnitude  and  lie  in  all 
quadrants  of  the  complex  plane;  therefore,  significant  numerical  round-off  error 
occurred  during  computation  of  the  residuals  within  algorithm  ACM  495 [  1  ] .  Even  if 
the  coefficients  of  best  approximation  had  happened  to  be  better  behaved,  serious 
cancellation  error  may  still  occur  in  some  problems  because  of  the  very  nature  of 
complex  arithmetic.  It  might,  therefore,  be  wise  to  use  a  double  precision  version  of 
algorithm  ACM  495  routinely  in  complex  Chebyshev  approximation  problems  io 
alleviate  such  cancellation  errors. 

One  method  of  detecting  the  presence  of  significant  round-off  errors  is  supplied 
by  the  nature  of  the  approximation  problem  itself.  That  is, Theorem  1  and  the  third 
step  of  the  proof  of  Theorem  2  together  imply  that 


Mnp(f)  ^  ^  Mnp(f)  SeC  (n/2P>  •  1 ) 

Once  Mnp(f)  and  the  coefficients  have  been  computed  in  algorithm  ACM  495,  these 
bounds  may  be  checked  to  see  if  significant  numerical  round-off  error  has  occurred. 
In  the  example  presented  in  the  Appendix,  rounding  to  5  significant  digits  gives 

.014436  =  Mnp(f)  <  d^p(f)  =  Mnp(f)  sec  (n/2p)  =  .014946  . 

However,  if  we  round  to  6  significant  digits  instead,  it  is  seen  that  the  second 
)  inequality  in  (6.1)  does  not  quite  hold.  We  conclude  that  the  effects  of  round-off 

errors,  although  visible  in  the  results,  are  not  significant  in  this  example.  (Single 
precision  numbers  on  the  DEC  VAX  11/780  have  approximately  7  significant 
decimal  digits.) 

A  sensitivity  analysis  on  the  optimum  coefficients  may  be  in  order  in  some  ap¬ 
plications  to  determine  their  utility.  This  consideration  is  completely  independent  of 
their  numerical  accuracy.  For  example,  in  an  antenna  array  design  problem  where 
some  elements  are  spaced  significantly  less  than  a  half-wavelength  apart,  it  might 
well  turn  out  that  the  optimum  coefficients  need  to  be  specified  with  a  relative  error 
of  better  than  10~6.  Then,  although  the  mathematical  results  may  be  correct  and 
accurate,  practical  usage  is  precluded.  This  sensitivity  can  be  determined  by  per¬ 
turbing  the  optimum  weights  a  few  percent  and  observing  if  a  drastic  change  occurs 
on  the  desired  sidelobe  behavior.  (Such  arrays  are  referred  to  as  super-directive 
arrays.) 
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Appendix 

The  Computer  Program 


$  RUN  ZACM 


CPU  TI*£  IN  SEC  s  1.61 
PAGE  FAULTS  s  1 
NURbEK  OF  ITERATIONS  *  17 


NUMBER  OF  DATA  POINTS  M 
NUMBER  OF  BASIS  FUNCTIONS  n 
NUMbER  OF  PHASES  PEk  POINT  P 


1^1 

3 

6 


LOWER  BGUnO  FOR  BEST  CHEBISrtEV  EPR"R 
UPPER  BOUND  FOR  BEST  CHe,B|Sri£V  EPAOk 
LOWER  BOUND  *  SECANT!  Pl/(2»P)  ) 


0.1 44  3&293E-0 1 
0.149460100-01 
0.i4945SSl£-yl 


CALCULATED  rank  *  6  , AS  EXPECTED. 

UNIQUE  SOLUTION: 


1 

2 

3 


REAL  PART 
0.37653J  93 
-2.012094/4 
2.64613126 


I.HAO  PART 
0.904025o7 
-2 . 01 409549 
1,09946144 


r>  o  no  on  on  o  o  orirjoooooooooonnnonooooonooootionnnooon 


t  TYPE  ZMAIN.FOH 


THIS  MAIN  ROUTINE  SOLVES  A  LINEAR  COMPLEX  FUNCTION 
APPROXIMATION  problem.  the  computed  approximation  can 
BE  MADE  AS  CLOSE  TO  THE  BtST  CHEPYSHEV,  OR  MlMMAX, 
APPROXIMATION  AS  DESIRED.  TnF  APPROXIMATION  IS  CONSTRUCTED 
On  A  FINITE  DATA  SET  FROM  ARBITRARY  BASIS  FUNCTIONS. 


REFERENCE: 

HOY  L.  STREIT  AND  ALBERT  «.  NUTTaLL,  "LInEAR  CHEBYSHEV 
COMPLEX  FUNCTION  APPROXIMATION,"  NUSC  TECHNICAL  REPORI  6403, 
NAVAL  UnOERmATER  SYST£VS  CEl.TER,  NE»i  LONDON,  CT, 06320. 


THIS  APPROACH  SOLVES  AT  MOST  M*P  LINEAR  EQUATIONS  In  2*N 
UNKNOWNS  IN  TnE  CrtEBYSHEV  NORM;  THAT  IS,  THE  MAXIMUM  MAGNITUDE 
RESIDUAL  IS  MINIMIZED.  ALL  EOUATIONS  ANU  UNKNOWNS  ARE  REAL. 

TnE  SOLUTION  IS  CQMPU'TFO  U5INU  LINEAR  PROGRAMMING. 

SINGLE  PRECISION  IS  USED  TO  SOLVE  THE  SYSTEM  OF  EQUATIONS; 
HOWEVER,  DOUBLE  PRECISION  IS  USED  TO  SFT  UP  THE  SYSTEM 
ITSELF  IN  ORDER  TO  MINIMIZE  POSSIBLE  ROUND-OFF  ERRORS. 
««««***•««***«*«**«********«*•* ******************************* 


THE  N  COMPLEX  COEFFICIENTS  OF  APPROXIMATION  ARE  GIVEN  BY* 
COEFOO  ♦  I  *  COFF(NfK)  ,  K»1,2,,..,N. 

ThE  COEF  ARRAY  IS  COMPUTED  IN  SUBROUTINE  ACM495. 


USERS  “UST  SPECIFY  THE  FOU-CWIKG  SIX  NUMBERS! 

ThE  NUMBER  OF  BASIS  FunCTIOwS: 

PARAMETER  N«3 

ThE  number  of  data  points: 

PARAMETER  “«101 

ThE  NuMbER  OF  PHASES  PFP  lA’CA  PUI..T:  [P  ,GF.  2] 

PARAMETER  P*6 

MUST  THE  FINAL  COEFFICIENTS  BE  REAL?  (iKEALsl  IFF  YES] 
PARAMETER  TREALaO 

MUST  THE  RESIDUALS  BF  MON-NEOATiVE?  ClSIOESal  IFF  YES] 
PARAMETER  ISIOES*0 

RELATIVE  ERROR  CRITERlO:  CmELFRRsO.O  IFF  CHEBYSHEV  SOLUTION] 
DATA  REuFRR/0,0/ 


r»  o  o  n  o  o  o  o  n  r>  o  n  n  o  r»  o  n  r>  o  o  n  r»  r>  n  n  n  o  r»  o 
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USERS  NEED  CHANGE  NOTh I NG  MORE  IN  THIS  ROUTINE,  BUT 
THEY  MUST  MAKE  APPROPRIATE  CHANGES  IN  THE  FOLLOWING 
SUBROUTINE'S  CALLED  SY  THIS  MAIN  ROUTINE! 


ZFUNCT 

ZSASIS 

ZABSCS 

ZPHASE 


DFFIi.ES  THE  FUNCTION  APPROXIMATED 
DEFINES  THE  BASIS  FUNCTIONS 
DEFINES  THE  FINITE  DATA  SET 
DEFIuES  THE  PHASES  AT  ALL  DATA  POINTS 


*******************************************¥***************** 


PARAMETER  MIPsM*P,MOIMsMIP*1,N2s2*N,NDTM*N2+3 
DIMENSION  BDATA(NDlM,MPIM),FDATAfMDIM) ,C0£F(NDIM)  , 

1  INDEX (  v) 

DOUBLE  PRECISION  COSDTA(MIP) , SInDT A CMlP) , ARGCMIP) ,ZRDATA(M) , 

1  ZIDATA(M),RESIDR(M) ,RESTOI(M) , CHEBER , DZERO, PI 

INTEGER  OCODE,RANK 
DATA  ZERO/O.O/ 

data  ozerq/o.oo/ 

DATA  PI / 1.1 41 5926535097 932 38DO/ 

SET  THE  UNIT  ROUND-OFF  ERROR  CFOR  VAX  11/780) 

DATA  TOL/.596E-7/ 

********** 

INITIALIZE 

********** 


IFCMIP.GE.NJG0  TO  50 
PRINT  51 ,MIP,N 

51  FORMAT ( / ,  '  *♦**♦  INITIALIZATION  ERROR !  M*P  MUST  EXCEED  N  '  , 

1  /,*  *****  aoTs  m*p  =  ',iio,',  and  n  *  ',110, 

2  /,'  ****»', 

3  /,*  *****  EXECUTION  TERMINATED') 

GO  TO  9999 
50  CONTINUE 

DO  89  1*1 ,N2 
COEF ( I ) =ZERO 
89  CONTINUE 

DO  91  1*1, M 
RtSIDR(I)* DZErO 
RESIDI (I )*DZEkO 
91  CONTINUE 

****************** 

DEFINE  THE  PROBLEM 
****************** 


DETERMINE  THE  DATA  POINTS 
CALL  ZASSCS(ZRDATA,ZIDATA,MJ 

DETlRMInE  THE  PHASES  AT  ALL  DaT«  po XNTS 
IP*P 

CALL  ZPHASFC  INDEX,  AUG,  ,  IP,  MlPSuM  ,  RtSIDR  ,  RESIDI) 
C 

C  COMPUTE  ThF  NECESSARY  SINES  AnD  COSINES 
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CALL  ZTRIGDt  ARG,  COSDTA,  STnDT A ,  MIPS’JM ,  INDEX) 

SET  UP  THE  OVER-DETERMINED  SYSTE*  OF  REAL  EQUATIONS 
CALL  ZFnSET (BDATA, FDAIA, COSDTA, Slh PTA, INDEX, NDIM,N,H, 

1  ZR DATA, ZIDAIA,  TREAD 

SET  CONSTRAINT  IF  COEFFICIENTS  MUST  BE  REAL 
NSETsn2 

IFdREAL.EQ.l  )N5EI*n 

SET  OPTION  FOR  ONt  SIDED  SOLUTION  OF  OVER -DETERM I NED  SYSTEM 
NSTDES*2 

IF ( IS1DES . EQ.  1 )NSID£Sal 

get  initial  timing  and  paging  information  (for  the  vax  h/7bo) 

CALL  GETJPKNCPUl  ,NPFS1) 

****************************************************** 

SOLVE  THF  OVER-DETERMINED  SYSYTEM  OF  mTPSUM  EQUATIONS 
TN  NSET  UnK»0aNS.  ALL  EQUATIONS  At.D  UNKNOWNS  ARE  REAL. 

******************************* **m* ********  *********** 

CALL  ACM495(MIPSUM,NS£TfM0lM,N0IM,BDATA,rDATA,T0L,RELERR,C0EF, 
1  RANK, RES* AX, ITER , OCODE , nSIOES ) 

COMPUTE  ThE  RESIDUALS  DIRECTLY  FOR  GREATER  ACCURACY 
CALL  ZRESID(R£SIOR,RFSIOI,N,M,COEF,ZRDATA,ZlDATA,CHEBER) 

GET  FINAL  TIMING  AND  PAGING  INFORMATION  (FOR  VAX  11/780) 

CALL  GETJPI(NCPU2,NPFS2) 

****************** 

PRINT  SUMMARY  DATA 
****************** 

PRINT  ELAPSED  TIMING  AND  PAGING  INFORMATION  (FOR  VAX  11/780) 
DCPU=(NCPU2-NCPU1)/100. 

IDPF=NPFS2-NPFSi 
PRINT  110,  PCPiJ 

110  FuRMAK/,'  CPU  TIME  IN  SEC  =  *,1F10.2) 

PRINT  U4,I0PF 

114  FORMAT! '  PAGE  FAULTS  =  *,IlO) 

PRINT  100, HER 

100  FORMAIC'  NUMBER  OF  ITERATIONS  =  ',110) 

PRINT  111,1 

111  FQPVAIC/,'  number  OF  DATA  POINTS  M  =  ',110) 

Pn InT  112, N 

112  Format c '  number  of  sasis  fuuCtiuns  n  a  ',iio> 

PRINT  115, P 

115  FORMAT ( '  NUMBER  OF  PHASES  PER  PJluT  P  a  ',110) 

IF  (IREAL.EQ.I)PRINT  134 

134  FORMAT!/,'  THE  COEFFICIENTS  ARF  REQoTRFD  TP  BE  REAL.') 

PRINT  103, RESMAX 

103  FORMA  I  ( / ,  '  LO.vER  SOuND  FOR  SEST  ChEoYSHEV  ERROR  a  *,1E16.8) 
PRInT  1  07  ,  CilEoER 

107  FORMAIC  UPPER  bOuNt)  f Oft  6E5T  PnFbYSuE V  ERROR  *  ',1016.8) 
SfcCaR£SMAX/COs(PI/(2.DO*P)) 
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PRINT  109, sec 

109  FORMAT!'  LOWER  BOUND  *  SEC*f.T(  PI/(2*P)  1  «  *,lElb.*) 

IF  (RANK.EQ.NSEnGO  TO  119 
PRINT  171, RANK 

171  FORM AT( / ,  '  CALCULATED  RANK  s  ',110,'  .') 

PRInT  118, NSFT 

118  FORMAT!/, S3! ***) ,/,53( ***),/, *  THE  kAf.K  SHOULD  EQUAL  ',110,*.', 

1  /,'  CHECK  FOR  POSSIBLE  ERRORS  IN  PROGRAM  ANu/OR  PROBLEM,', 

2  /,5J(**'),/,53('*')) 

GO  TO  181 

119  CONTINUE 
PRINT  117, RANK 

117  FORMAT!/,*  CALCULATED  RANK  *  *,110,*  ,  AS  EXPECTED.*) 

181  CONTINUE 

PRINT  RELATIVE  ERkOn  OF  MAXI*!)*  RESIDUAL  IF  NECESSARY 
IF(KELERR.GT.Z£RO)PkInT  124, RELERR 
124  FORMAT!//,*  RELATIVE  ERROk  IN  THE  MAXI*U-  RESIDUAL  INCURRED*, 

1  /,*  BY  THIS  APPROXIMATE  SOLUTION  *  *,1£16.8,/) 

IF(OCODE.EQ.O)PRlNT  121 
IFI0C0DE.EQ.1 )PRTnT  122 
IF!0C0D£.EQ.2)PkI,nT  123 

121  FORMAT!/,'  THE  FOLLOWING  SOLUTION  IS  PROMAPLY  NON-UNIQUE »*,/ ) 

122  FORMAI ! / , *  UNIQUE  SOLUTIONS',/) 

123  FORMAI!/,'  PREMA'f JR£  TERMINATION  DUt  10  ROUND-OFF  ERRORS.', 

I  /,'  BEST  COMPUTED  SOLUTION:',/) 

PRINT  133 

133  FORMATI15X,  *R£Al,  PART',  7x,'ImAG  PAkT') 

*************************************** 

PRINT  THE  COEFFICIENTS  OF  APPROXIMATION 
*************************************** 

PRINT  102,  !I,COEF(I),cnEF!M+I),l*i,,,) 

102  F0RMAT!I5,3X,2F16.8) 

******************* 

EnD  OF  main  PROGRA-' 

******************* 

9999  END 


TYPE  zfunct.for 


SUPROyTlNt  E* UVCT  EVALUATES  Tn F  APPROXIMATED 
CuMPL£X  FUNCTION  AT  A  SPECIFI&D  Data  POIVI. 

********************************************* 

ZFNR  s  REAL  PART  uF  TmF  FuiCTiO*  (OUTPUT) 

ZFNi  s  1MAG  PART  OF  THt  FuNCTIO*  (OuTPUT) 
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ZR  =  REAL  PART  OF  TnE  DATA  POINT  (INPUT) 
ZI  s  I*AG  PART  OF  THE  DATA  POINT  (INPUT) 

tv******************************************* 

SUBROUTINE  ZFUNCT( Zf NR , ZFNI , ZR , ZI ) 

DOUBLE  PRECISION  ZR , ZI , ZFhR , Zr N1 

************************************ 

BEGIN  USFR  CODE  FOR  SPECIFIC  PROBLEM 
************************************ 

ZFNR*C0S(3.D0*ZR) 

ZFNI*STN(3.00*ZR) 

********************************** 

END  USER  CODE  FOR  SPECIFIC  PROBLEM 
********************************** 


************************ 
END  OF  SUBROUTINE  ZFUnCT 
************************ 

RETURN 

END 


TYPE  ZBASIS.FOR 


********************************************* 

SUBROUTINE  ZBASIS  EVALUATES  THE  IB-TH  COMPLEX 

T. 

****** 

************** 

(INPUT) 

CTIO\,  (OUTPUT) 
ctiom  (Output) 
POINT  (I.VPuT) 
POINf  (INPUT) 

************** 

SUBROUTINE  ZRmSIS i I e , ZF  «B , ZF"I »  ZR  » 21 ) 

Double  precise  zR,zi,zFt.s,zr*<’i 

************************************ 

«EGIN  tloER  COLE  Fur  SpFCIrTC  pR0°UEB 


BASIS  FUNCTION  AT  a  SPECIFIED  DATA  POIN 

*************************************** 


THE  BAjIS  FUNCTION  I„DtX,  *»He.P£ 
IB  *  t,  l.  ... 

REAL  PART  QF  ▼  „£  T a-Ty  BASIS  BUN 
IMaG  PART  OF  T„F  lo-TH  BASIS  FWN 
REAL  PART  OF  TriF  SPECIFIED  DA?A 
I**G  PART  OF  TnE  SPECIFIED  Data 


7F*ik=CC!S(  (  Td-i  ,00)*zR; 
ZFNl*SlN((I6-l.uO)*ZR) 
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END  USER  COOE  FOR  SPECTFIC  PROBLEM 

********************************** 


************************ 
END  OF  SUBROUTINE  Z8ASIS 
************************ 

RETURN 

END 


TYPE  ZAbSCS , FOR 


******************************************** 
SUBROUTINE  ZABSCS  DEFINES  THE  DATA  POINTS  ON 
WHICH  THE  APPROXIMATION  IS  CONSTRUCTED. 

******************************************** 

***************************************************** 

ZRDATA(I)  a  REAL  PART  OF  THE  I-IH  DAtA  POINT  (OUTPUT) 

ZIDATACI)  a  ImAG  PART  CF  THE  I-TH  DATA  POINT  (OUTPUT) 

N  a  TOTAL  NUMBER  OF  DaTA  POINTS  (INPUT) 

***************************************************** 

SUBROUTINE  ZABSCS(ZRDATA,ZIDATA , M ) 

DOUBLE  PRECISION  ZRDAT A ( l ) , ZIDAIA (  1 ) 

************************************ 

BEGIN  USER  CODE  FOR  SPFCTFIC  PROBLEM 
************************************ 

DOUBLE  PRECISION  X,PI 

DATA  PI/3.U15926S35R9793238D0/ 

DO  10  Isl,M 

X=PI*(T-1  .D0)/C4.D0*(M-1.D(>)  ) 

ZRDATA ( I ) sX 
ZIDATA ( I ) aO . DO 
10  CONTINUE 

********************************** 

End  USER  CODE  FOR  SPECIFIC  PROBLEM 

********************************** 


************************ 
End  OF  SUaPuUIIi.E  ZABSCS 
************************ 

return 

END 
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$  TYPE  ZPHASE.FOP 


SUBROUTINE  ZPnASE  DEFINES  ThE  PriASES  AT  FACH  DATA 
POINT,  AS  KELL  AS  THE  TOTAL  NUMBER  OF  PHASES. 

NOTE  THAT  ’’’riE  TOTAL  NUMBER  OF  PHASES  EOuALS  THE 
TOTAL  NUMBER  OF  EuI'aTIOnS  SOLVED  BY  ACM495. 

************************************************* 

***  CAUTION  ***  THE  TOTAL  NUMBER  OF  PHASES  MUST  NEVER 
***************  EXCEED  the  product  M*P.  however,  the 
***************  NUMBER  OF  PHASES  MAY  BE  VARIED  FROM  DATA 
***************  POINT  TO  OATA  POINT.  SEE  NUSC  REPORT  6043. 

*************9***9************«******«******«**************** 

INDEX ( I )  »  NUMBER  OF  PHASES  AT  THE  I-TH  DATA  POINT  (OUTPUT) 

ARG  a  AGGREGATE  ARRAY  OF  ALL  PHASES  AT  ALL  DATA 

POINTS  IN  LEXICOGRAPHICAL  ORDER  (OUTPUT) 

M  a  NUMBER  OF  OATA  POINTS  (InPUT) 

P  a  THE  INTEGER  PARAMETER  P  OF  MAIN  ROUTINE  (INPUT) 

MIPSUM  a  TOTAL  NUMBER  OF  PHASES  IN  THE  ARRAY  ARG  (OUTPUT) 

THE  FOLLOWING  TWO  ARRAYS  ARE  ZFRO  FILLED,  UNLESS  AN 
INITIAL  APPROXIMATION  HAS  BEEN  PROVIDED  IN  THE  MAIN 
ROUTINE.  USE  ONLY  IF  NFEDED,  OR  ELSE  IGNORE  THEM. 

RESIDR ( I ) a  HEAL  PART  OF  CuMPLEX  RESIDUAL  AT  THE  I-TH 

DATA  POINT  (INPUT) 

FESIOI(I)a  IMAG  PART  OF  COMPLEX  PESIDUAl  AT  THE  I-TH 

OATA  POINT  (INPUT) 

************************************************************ 

SUBROUTINE  ZPHASE ( INDEX , ARG, H , P ,« IPSUM , RESIDR , RESIDI ) 

integer  p 

DIMENSION  I n»D£X (  1  ) 

DOUBLE  PRECISION  ARG ( 1) , PESIDK f I) , PESI DI ( 1 ) 

************************************ 

BEGIN  USER  CODE  FOP  SPECIFIC  rRu»LE« 
************************************ 

THIS  CODE  DEFINES  TrtF  PhASES 

pi*(j-n/p  ,J*t,?,...,p, 

AI  EACH  OF  iHE  m  DATA  POINTS. 

DOUBLE  PRECISIOi.  PI,X 
DATA  PI/3.l4l5926‘>35»979323»D0/ 

DEFINE  THE  NUMBER  Of  PHAStS  AT  EACH  DATA  POINT 

DC  1"  I a  1  ,  M 
lNOEX(I)ap 
10  CONTINUE 
C 

C  OfcFINE  ALL  THE  PHASES 


34 


noonn  o  o  o  o  n  n  m  nnnnnnnnnn  nnn 


ML00P*0 
DO  20  1*1, M 
LOOP* INDEX C I ) 

X*L00P 

DO  30  J*1 , LOOP 
ML00P*ML00P+1 
ARG(MLOQP)*PI*( J*l)/X 
30  CONTINUE 
20  CONTINUE 

TOTAL  NUMBER  OP  PHASES 

MIPSUM*ML00P 

***************9****************** 

END  USER  CODE  FOR  SPECIFIC  PROBLEM 

********************************** 


************************ 
END  OF  SUBROUTINE  ZPHASE 
************************ 

RETURN 

END 


TYPE  ZTRIGD. FOR 


********************************************************* 
SUBROUTINE  ZTRIGD  COMPUTES  THE  REQUIRED  SINES  AND  COSINES 
*********************** a********************************* 

SUBROUTINE  ZTRIGDURG,  COSDTA,  5INDTA,MIPSUM, INDEX) 

DOUBLE  PRECISION  X, ARG (1 ) , COSDTA ( 1 ) , SINDTA ( 1) 

DIMENSION  INDEX  Cl) 

DO  10  1*1, M IPSUM 
XsARG(I) 

COSDTA ( I)sCOS(X) 

SINDTA(I)*SIN(X) 

10  CONTINUE 

************************ 

END  OF  SUBROUTINE  ZTRIGD 
************************ 

RETURN 

END 


«  TYPE  ZFN5ET . FOR 

i 
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I********************************************************* 

SUBROUTINE  ZFnSET  SETS  UP  THE  TRANSPOSE  nF  THE  COEFFICIENT 
MATRIX  OF  THE  REQUIRED  OVER-DETERMINED  SYSTEM  OF  EQUATIONS 
IN  THE  ARRAY  NAMED  riDAT A .  SPECIFICALLY, 

TRANSPOSECROATA)  *  a  =  FDATA 

WHERE  FDATA  IS  THE  CONSTANT  VECTOR  AND  X  IS  THE  UNKNOWN 
SOLUTION  VECTOR.  THE  FDATA  ARRAY  IS  ALSO  FILLED  BY  THIS 
SUBROUTINE.  IF  THE  X  VECTOR  IS  RFQUIRED  TO  BE  REAL,  THEN 
ThE  "LAST"  HALF  OF  EACH  Or  THE  REQUIRED  EQUATIONS  ARE 
IGNORED  IN  THE  FINAL  COMPUTATIONS.  HENCE,  THIS  PART  OF  THE 
BDATA  ARRAY  IS  NOT  COMPUTED  IF  AND  ONLY  TF  IREAL  =  1  . 

**«V*«***«*«*f*«*****f*************************»******«*** 

SUBROUTINE  ZFNSET ( BuAT A , FD AT A , CUSDTA , SI  NOT A , INDEX , NDlM , N , M , 

1  ZRDATA,ZIOArA,IREAL) 

DIMENSION  INDEX  (  1)  ,  BOAT  AUDI*,  1)  ,Cf)ATA(  1 ) 

DOUBLE  PRECISION  ZH , Zi , ZRO AT A ( 1 ) , ZI 0 A1 A C 1 ) , COSDTA ( 1 ) ,SINDTA(1) , 
1  ZFwR , ZFN I 

FILL  THE  FDATA  ARRAY 

MLOOPsO 
DO  10  1  =  1,  M 
ZRsZRDATAClJ 
ZI  =  ZI 0  AT  A ( I ) 

CALL  ZFUNCT(ZFNR,ZEnI,ZR,4I) 

LOO?=INOEX( I) 

DO  20  J=1 .LOOP 

mloop=mloop*i 

FDA TA(mlOOPJ=ZFnR*COSDTA(MLU0P)+ZFm 1*51 nDIA(MLOOP) 

20  CONTINUE 
10  CONTINUE 

FILL  THE  BDATA  ARRAY 

00  JO  K=1,N 
IB  =  K 
WLOOP=0 
Du  40  1  =  1,  m 
ZRsZRDATA(I) 

ZI=ZIDATA(I) 

CALL  Z«ASIS(IS,ZF,'.R,ZFNl,ZR,Zn 

LGO?=INOEX(n 

Du  SO  J=1 , LOOP 

MuOuP=MuOQ°+t 

BuAIACK,MUOOB)=ZF\R*COSDTACMlOOp)+Zr N 1  * Si NOT A C MLOOP) 

TF  (TREAL.FQ. I )G0  TO  50 

BDAIAIK  +  N,  MbOOP)=ZFNR*SINLTAC,VLOnP)-ZfNI*COSDTACMLOOP) 

SO  CONTINUE 
40  CONTINUE 
JO  CONTINUE 


c 
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END  OF  SUBROUTINE  ZFNSET 

************************ 

return 

END 


TYPE  ZRESID . FOR 


******************************************************* 
SUBROUTINE  ZRESID  COMPUTES  THE  COKPlEX  RESIDUAL  AT  EACH 
OF  THE  OATA  POINTS  DIRECTLY  TROM  THfc  DEFINITION. 

******************************************************* 

********************************************************** 

RESIDR C I )  =  REAL  PAkT  OF  RESIDUAL  AI  THE  I-TH  OATA  POINT 

RESIDICI)  =  IMAG  PART  OF  RESIDUAL  AT  THE  I-TH  DATA  POINT 

N  =  HUMBER  OF  BASIS  FuNCTIOkS 

M  a  NUMBER  OF  DATA  POINTS 

COEF(I)  *  REAL  PART  OF  THE  I-TH  COElFIClEfcT,  1*1,2, ...,N 

COEFCN  +  I)  *  IMAG  PART  OF  THE  1-1H  COEFFICIENT,  1*1,2,  ...,N 

ZRDATA ( I )  a  REAL  PART  OF  THE  I-TH  DATA  P<MNI  ,  1*1,2,...,* 

ZIOATACI)  =  I*AG  PART  OF  THE  I-TH  DATA  POINT  ,  I*1,2,...,M 

CHEBER  a  HAXIMu*  *AGnIIUCE  RfcSIDOAL  l rHEbYSHtV  ERROR) 
********************************************************** 

SUBROUTINE  ZRESID ( RESIDR , RES^Dl , K ,**, CCFF , ZRDATA , ZIDATA , 

1  CHESER) 

DOUdLE  PRECISION  RESlUR ( 1 ), RESIDIC 1 ), ZRDATA ( 1 ), ZIDATA ( 1 )  , 

1  ZR,Z  I,  ZFnR,ZFM,  CHEBER,  RESR.RESI 

DIMENSION  COEF(l) 

C 

CHEBERa-1 .DO 
DO  10  1*1, M 
ZRsZRDATA(I) 

ZI*ZIOAI A ( I ) 

C 

CALL  ZFLNCT(ZFNR,ZFNI,ZR,ZI) 

c 

RESR=ZFNR 
RESI=ZFN I 
ZRsZRQA IA ( I ) 

ZI*ZIOATA(I) 

DO  20  J*1,N 
IBSJ 
C 

CALo  ZBASlSlIB,ZFrtR,ZF«I,ZP,Zl> 

C 

CuEFRsCOEF ( U ) 

COEr IaCCEF (N+J) 

RESRaREoP-C COEFR*ZF nR-COFFT*ZF  '* 1 5 
R£Sl*REST-(COt'FR*ZFM  +  CuEFI*ZFNR) 

20  CONTINUE 

R£SIDR(I)*R£SR 
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RESIDICIJaRESI 

C 

ZR*SQRT(RESR*RSSR+RESi*RESI) 
IF(ZR.GT.CH£B£R)CHEBERaZR 
10  CONTINUE 


************************ 
End  oe  subroutine  zresid 
************************ 

RETURN 

END 


TYPE  ACM495.F0R 


************************************************************ 

SUBROUTINE  ACM495  SOLVES  THE  REAL  OVER-DETERMINED  SYSTEM  OF 
M  LINEAR  EQUATIONS  IN  N  UNKNOWNS 

TRANSPOSE(A)  *  X  *  fc 

IN  THE  CHEBYSHEV  NORM;  THAT  IS,  THE  MAXIMUM  MAGNITUDE 
COMPONENT  OF  THE  RESIDUAL  VECTOR  t  a  B-TRANSPOSE(A)*X  J 
IS  MINIMIZED  OVER  ALL  CHOICES  OF  THE  UNKNOWN  VECTOR  X. 

ONE  SIDED  CHEBYSHEV  SOLUTIONS  CAN  ALSO  BE  COMPUTED;  THAT  IS, 

THE  MAXIMUM  COMPONENT  OF  THE  RESIDUAL  VECTOR  IS  MINIMIZED 
SUBJECT  TO  THE  REQUIREMENT  THAT  ThE  RESIOUAL  VECTOR  BE 
NON-NEGATIVE,  CSEE  REF.  2, PAGE  863.] 

************************************************************ 

*************************************************************** 

REFERENCE  S 

1.  I.  BARRODALE  A HO  C.  PHILLIPS,  "SOLUTION  OF  AN  OVERDETERMINED 
SYSTEM  Of  LINEAR  EQUATIONS  IN  THE  CHEBYSHEV  NORM,"  ACM 
TRANS.  ON  MATH.  SOFTWARE,  VOL.  1,  NC.  3,  SEPTEMBER  1975, 

?P.  264-270. 

2.  I.  BARRODALE,  L.  M.  DELVES,  A NO  J.  C.  MASON,  "LINEAR 
CHEBYSHEV  APPROXIMATION  OF  Cu* PLEX-VALoED  FUNCTIONS," 

MATH.  CnMP .  ,  VOL.  32,  .*0.  143,  JULY  197S,  PP.  863-863. 

*************************************************************** 

SUBROUTINE  ACM49S C M , 4 , MDlM , ft DIM , A , B , TOL , OELERR , X , R ANK , 

1  RESMAX,ITcP,OCODE,NSIuES) 

*********************************************************** 

M  a  NUMBER  OF  EQUATIONS. 

N  =  NUMBER  OF  UNKnO*N5.  CN  MUST  NOT  EXCEED  M) 

“DIM  a  NUMBER  OF  COLUMNS  OF  A,  (MDIM.GE.M+1) 

NOIM  a  NUMBER  Of  ROwS  UF  A,  CNDIm ,GE . V+3 ) 

A  a  T*0  DIMENSIONAL  REAL  ARRAY  DlMFnSlONFD  A(NDIM,MDIM) . 

ON  ENTRY,  THE  TRANSPOSE  OF  THE  MATRIX  OF  COEFFICIENTS 
OF  THE  OVER-DETERMINED  SYSTE*  MuSI  BE  STORED  IN  THE 
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FIRST  V  COLUMNS  AND  N  ROWS  OF  A.  THESE  VALUES  ARE 
DESTROYED  BY  THIS  SUBROUTINE. 

8  *  ONE  DIMENSIONAL  REAL  ARRAY  OF  LENGTH  MDIM. 

ON  ENTRY,  8  CONTAINS  THE  RIGHT  HAND  SIDES  OF  THE 
M  EQUATIONS  OF  THE  OVER-DETERMINED  SYSTEM  IN  ITS  FIRST 
M  LOCATIONS.  ON  EXIT,  B  CONTAINS  THE  RESIDUALS  FOR  THE 
EQUATIONS  IN  ITS  FIRST  M  LOCATIONS.  (SEE  NSIDES.)  THESE 

residuals  are  not  computed  directly  from  the  definition. 

T0L  *  A  SMALL  POSITIVE  TOLERANCE,  TYPICALLY  THE  UNIT 
round-off  error  of  the  computer. 

RELERR  >  A  REAL  VARIABLE  WHICH  ON  ENTRY  MUST  EQUAL  0. 

IF  A  CHESYSHEV  SOLUTION  IS  REQUIRED.  IF  RELERR 
IS  POSITIVE,  THIS  SUBROUTINE  CALCULATES  AN 
APPROXIMATE  SOLUTION  WITH  RELERR  AS  AN  UPPER 
BOUND  ON  THE  RELATIVE  ERROR  ON  ITS  LARGEST 
RESIDUAL,  ON  EXIT,  RELERR  GIVES  A  SMALLER  UPPER 
BOUND  FOR  THIS  RELATIVE  ERROR.  [SEE  REF.  1.] 

X  s  ONE  dimensional  REAL  ARRAY  OF  length  NDIM. 

ON  EXIT,  X  CONTAINS  A  SOLUTION  TO  THE  PROBLEM  IN  THE 
FIRST  N  LOCATIONS.  THE  CONTENTS  OF  X(N*i) , . . . ,X(NDIM) 

ARE  UNCHANGED. 

RANK  *  AN  INTEGER  WHICH  GIVES  ON  EXIT  THE  RANK  or  THE 
COEFFICIENT  MATRIX.  (WILL  DEPEND  ON  TOL.) 

RESMAX  b  ON  EXIT,  THE  LARGEST  MAGNITUDE  RESIDUAL. 

OCODE  b  0  IF  OPTIMAL  SOLUTION  IS  PROBABLY  NONUNIQUE.  A  DEFINITE 
STATEMENT  REQUIRES  FURTHER  COMPUTATION  WHICH  IS  NOT 
DEEMED  to  BE  COST  EFFECTIVE. 

»  1  IK  UNIQUE  OPTIMAL  SOLUTION. 

a  2  IK  CALCULATIONS  TERMINATED  PREMATURELY  DUE  TO 
ACCUMULATED  ROUND-OFF  ERRORS. 

NSIDES  s  1  IF  ONE  SIDED  CHESYSHEV  SOLUTION  IS  COMPUTED. 

IN  THIS  CASE  THE  RESIDUALS  RETURNED  FROM  THIS 
PROGRAM  ARE  ERRONEOUS  AND  MUST  BE  COMPUTED  IN 
THE  CALLING  ROUTINE.  [SEE  REF.  2.1 
s  2  IF  TWO  SIDED  CHEBYSHEV  APPROXIMATION  IS  COMPUTED. 

THIS  IS  THE  STANDARD  FORM,  THE  RESIDUALS  RETURNED 
FROM  THIS  PROGRAM  ARE  CORRECT  IN  THIS  CASE.  MORE 
NUMERICAL  ACCURACY  IN  THE  RESIDUALS  MAY  RESULT 
FROM  DIRECT  CALCULATION  IN  THE  CALLING  PROGRAM. 
******* ********************************************************* 

DIMENSION  A (NDIM, MDIM) ,B(MDIM) ,X (NDIM) 

INTEGER  PROW, PCOL, RANK, RANKP1 , OCODE 

THE  FOLLOWING  NUMBER  IS  MACHINE  DEPENDENT. 

DATA  BIG/1. £+3«/ 

************** 

INITIALIZATION 

************** 

IF (N SIDES. LE. 1 )SIDES=1 . 

IF(nSIDES.GE.2)5IDESb2. 

MPIbM+1 

NPlsN+1 

NP2sN*2 

NP3»N*3 

NPlMRxl 

RANK»N 
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RELTMP*RELERR 

RELERR*0. 

DO  10  J«1,H 
A(NP1,J)«1. 

A(NP2,J)*-B(J) 

A(NP3,J)*N+J 
10  CONTINUE 

A(NP1,MP1)»0. 

ITER*0 
OCODEsl 
DO  20  1*1 , N 
X(I)«0. 

A(I,«Pl)*I 

20  CONTINUE 

******* 

LEVEL  1 
******* 

LEV*1 

K*0 

30  K*K*1 

THE  NEXT  VARIABLE  IS  NOT  USED,  BUT  IS  IN  THE  PUBLISHED  CODE. 
KP1*K*1 

NPIMKsNPl-K 
MOD£*0 
DO  40  J=K,M 
B(J)*1. 

40  CONTINUE 

DETERMINE  THE  VECTOR  TO  ENTER  THE  BASIS 

SO  D=-BIG 

DO  60  JsK,M 

IF(BCJ) .EQ.O.JGO  TO  60 
DDsABS (A(NP2,J)) 

IF(DD.LE.D)GO  TO  60 

PCOLsJ 

D=DD 

60  CONTINUE 

IFU.GT.UGO  TO  70 

TEST  FOR  ZERO  RIGHT  HANO  SIuE 

TFCD.GT.TOLJGO  TO  70 
RESMAXsO. 

MOD£»2 
GO  TO  380 

DETERMINE  THE  VECTOR  TO  LEAVE  ThE  *ASI5 

70  DsTOL 

DO  80  I*1,N?1MK 
DD*A8S  C  A ( I , PCOL) ) 

If (OD.LE.07GO  TO  60 
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80  CONTINUE 

IFCD.GT.TOLJGO  TO  330 

CHECK  FOR  LINEAR  DEPENDENCE  IN  LEVEL  1 
B(PC0L)»0. 

IFCM00E.EQ.DG0  TO  80 
.  00  100  J*K,M 
IFCB(J).EQ.0.)GO  TO  100 
DO  90  I«1,NP1MK 

IF ( ABS ( A( I , J ) ) t LE.TOL) GO  TO  90 
MODEsl 
GO  TO  50 
90  CONTINUE 
100  CONTINUE 
RANK*K-1 
NPIMRaNPl-RANK 
OCODEsO 
GO  TO  160 

110  IF(PC0L.EQ.K)G0  TO  130 

INTERCHANGE  COLUMNS  IN  LEVEL  1 

DO  120  1*1, NP3 
D=A(I , PCOL) 

A  ( I , PCOL  )*A(I,K) 

A(I,K)*D 
120  CONTINUE 

130  IFfPROW.EQ.NPlMKJGO  TO  150 

INTERCHANGE  ROWS  IN  LEVEL  1 

DO  140  J«1,MP1 
D*A(PRON,U) 

A(PR0», J)*A(NP1MK, J) 

A(NP1MK,J)*0 
140  CONTINUE 
150  IFCK.LT.NJOQ  TO  30 

160  IFCRANK.EQ.iOGO  TO  380 

RANKP1*RANK+1 

******* 

LEVEL  2 
******* 

L£V*2 

DETERMINE  THE  VECTOR  TO  ENTEP  THE  BASIS 
D*T0L 

DO  170  J*RA«KP1,m 
PD*ABSCACaP2,J)) 

IF ( DD , L£ , D ) GO  TU  170 

PC0L*U 

D*DD 

170  CONTINUE 

COMPARE  CHEBYSHEV  ERROR  WITH  TOL 
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IF(D.GT.T0L)G0  TO  180 
RES MAX*0. 

M0D£*3 
GO  TO  380 

180  IF(A(NP2,PC0L).LT.-T0L)G0  TO  200 
A ( NP1 ,  POOL )> SIDES* A (NP1 , PCOL) 

DO  190  I*NP1 MR , NP3 
IFCI .E0.NP1 )G0  TO  190 
A(I,PCOL)»-A(I,PCOL) 

190  CONTINUE 

ARRANGE  FOR  ALL  ENTRIES  IN  PIVOT  COLUMN 
(EXCEPT  PIVOT)  TO  BE  NEGATIVE 

200  DO  220  IsNPIMR.N 

IF(A(I,PCOL),LT.TOL)GO  TO  220 
DO  210  Jsl,M 

A(NP1 , J)«ACNP1,J )f SIDES* Ad,  J) 
ACI,J)*-A(I,J) 

210  CONTINUE 

A(I,MPl)a»A(I,MPl) 

220  CONTINUE 
PROMsNPl 
GO  TO  330 

230  IF(RAhKPl.EQ.M)GO  TO  390 
IF ( PCOL . EQ . M ) GO  TO  250 

INTERCHANGE  COLUMNS  IN  LEVFL  2 

DO  240  I«NP1MR,NP3 
D*A (I » PCOL) 

Ad  ,PCOL)*A(I,M) 

A (I #  M )»D 
240  CONTINUE 
250  MM1«M-1 

******* 

LEVEL  3 
******* 

LEV  =  3 

DETERMINE  The  VECICR  TO  ENTER  ThF  BASIS 

260  Ds-TOL 

VAL=SIDES*A(NP2,M) 

DO  280  JsRAi.KPl  ,MMl 
IFCA(NP2,J) ,5ft. D)GC  TO  270 
PCOLsj 
D*A(NP2,J) 

M0DE*0 
GO  10  280 

270  DD»VAL-A(NP2, J) 

IF (DD.GE.D)GO  TO  280 

MGDE»1 

PCOL*J 

D*DD 

280  CONIINUE 

IF(O.GE.-TOL)GO  TO  380 
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DD*-D/A(NP2,M) 

IF(DD.GE,RELTMP)GQ  TO  290 

RELERMDD 

M0D£*4 

GO  TO  380 

290  IF(MODE.EQ,0)GO  TO  310 
DO  300  I*NP1MR,NP1 
A C I, POOL) »S IDES* A ( I, M) -A (I, PC0L) 

300  CONTINUE 

A(NP2,PCOL)*D 
A ( NP3 , PCOL) **A(NP3 , PCOL) 

DETERMINE  THE  VECTOR  TO  LEAVE  THE  BASIS 

310  D=BIG 

DO  320  I*NP1MR#NP1 

IF ( A ( I , PCOL} ,LE,TOL)GO  TO  320 

D0»A(I,M)/A(I,PC0L) 

IF C DD.GE . D) GO  TO  320 

PROw*I 

D*DD 

320  CONTINUE 

IF(D.LT.BIG)GO  TO  330 

0C0DE*2 

GO  TO  380 

PIVOT  ON  A(PRQW,PCOL) 

330  PIVOT»A ( PROW , PCOL) 

DO  340  J«1,M 

A (PROM, J)«ACPR0*,J) /PIVOT 
340  CONTINUE 

DO  360  J«1,M 

IF ( J « £0 . PCOL ) GO  TO  360 

OsA(PRON,J) 

DO  350  I«NP1MR,NP2 
IF(I.£Q.PRON)GO  TO  350 
ACI,J)sA(I,J)-D*ACI,PCOL) 

350  CONTINUE 
360  CONTINUE 

TPIVOTs-PIVOT 
DO  370  I»NP1*R,NP2 
A  ( I , PCOL) * A  C 1 , PCOL) /TPI VOT 
370  CONTINUE 

A (PRO*, PCOL) »1. /PIVOT 
D»A (PRO* , MP1 ) 

A(PR0rf,MPl)sA(NP3,PC0L) 

A(NP3,PCOL)aO 

ITER*IT£R+1 

GO  TO  (110, 230, 260), LEV 

************** 

PREPARE  OUTPUT 
************** 

380  DO  390  J«l,rt 
8 ( J ) *0 . 

390  CONTINUE 

If (M00E.EQ,2)G0  TO  450 


nonod  <"»  n  o 


TR  6403 


DO  400  0*1, RANK 
K*ACNP3,J3 
XCK)*ACNP2, J) 

400  CONTINUE 

IFCMODE.£Q.3.QR.RANK.EQ.>OGO  TO  450 
DO  410  I*NP1MR,NP1 
K*ABS(A(I,MPm-PLOAT(N) 
BCK)*ACNP2,M)*SIGNC1.,ACI,MP1)) 

410  CONTINUE 

IFCRANKP1.EQ.JOGO  TO  430 
DO  420  J*RANKP1,NM1 
KsABS(A(NP3,0) ) -FLOAT CN) 

BCK)*(ACNP2,M)-ACNP2,jn*SIGNCl .  ,ACNP3,J1) 
420  CONTINUE 

TEST  FOR  NON-UNIQUE  SOLUTION 

430  DO  440  I*NP1NR,NP1 

If CABSCACI,M>) ,GT.TOL)GO  TO  440 
OCODEsO 
GO  TO  450 
440  CONTINUE 

450  If (rtOOE.NE.2.ANO.MOOE.NE.33RESNAX*A(NP2,M) 
If (RANK. EQ.MlRESMAXsO. 

IF CB0DE.EQ.4)RESMAX*RESKAX-0 

************************ 

END  OF  SUBROUTINE  AC»495 
************************ 

RETURN 

END 
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